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A  multi-block  computational  method  is  developed  for  three-dimensional 

high  speed  turbulent  flows  over  complex  configurations.  In  this  method,  the  flow 

domain  is  divided  into  several  contiguous  blocks  such  that  each  block  is  partially 

bounded  by  a  soUd  surface.  The  chosen  grid  topology  has  advantages  in 

computing  eddy  viscosity  as  well  as  in  applying  a  thin-layer  Navier-Stokes  code. 

In  the  solution  process,  each  block  is  then  treated  as  an  independent  flow 

problem  governed  by  the  same  set  of  time  dependent  thin-layer  Navier-Stokes 

equations.  The  interface  boundary  conditions  between  blocks  are  updated  at 

every  time  step  of  the  solution  algorithm.  The  application  of  the  method  with 

distance  weighted  interpolation  for  updating  the  interface  boundary  condition  is 

rather  simple  and  straightforward;  however,  special  measures  in  grid  topology  and 


V 


generation  are  required.  The  flow  solver  employed  is  a  total  variation 
diminishing  thin-layer  Navier-Stokes  code  implemented  with  an  algebraic  eddy 
viscosity  model  of  Baldwin-Lomax.  In  this  study,  the  proposed  method  is 
investigated  for  the  simulation  of  a  transonic  turbulent  flow  past  a  wing-fuselage 
configuration  in  which  the  flow  field  is  properly  divided  into  six  blocks.  A  coarser 
grid  was  first  used  for  flow  simulation,  followed  by  grid  refinement  smdies.  With 
the  use  of  the  refined  grid,  the  calculated  wing  surface  pressure  distribution 
agrees  well  with  measured  data.  Numerical  results  obtained  also  show  that  the 
computed  shock  location  and  pressure  distribution  on  the  wing  can  be  in  good 
agreement  with  the  experiment  data  if  the  block  grids  used  are  adaptive  to 
important  solution  characteristics.  It  is  concluded  that  the  proposed  method  is 
indeed  a  very  promising  approach  to  be  developed  further  for  complex 
configuration  flow  simulation. 
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CHAPTER  I 
INTRODUCTION 


For  a  high  speed  flight  vehicle,  numerous  complex  three-dimensional  fluid 
flow  phenomena  appear,  including  the  formation  and  shedding  of  vortices,  the 
shock-wave  boundary-layer  interaction  and  induced  flow  separation,  and  the  shock 
on  shock  interaction.  Since  the  performance  of  a  flight  vehicle  is  strongly 
affected  by  these  complex  flow  phenomena,  accurate  prediction  of  them  and 
associated  aerodynamic  forces  and  moments  are  critical  for  advanced  design  of  a 
modem  flight  vehicle.  At  present  an  accurate  wind-tunnel  simulation  of  the 
complex  flow  problem  in  flight  conditions  can  be  very  difficult,  if  not  impossible. 
With  the  advent  of  supercomputers  and  the  advancement  of  solution  techniques 
for  nonlinear  problems,  computational  fluid  dynamics  (CFD)  is  under  extensive 
development  and  has  been  apphed  to  complement  the  wind-tunnel  experiment 
and  field  test  in  the  design  process  of  a  flight  vehicle. 

The  use  of  CFD  in  the  aerodynamic  design  of  flight  vehicles  has  / 
progressed  rapidly  over  the  last  few  years.  In  the  early  days,  CFD  was  used  to 
support  the  vaUdity  of  a  design  that  developed  in  the  wind-tunnel  by  trial  and 
error.  The  state-of-the-art  of  the  CFD  method  has  progressed  to  the  point  of 
being  regarded  as  an  important  design  tool  for  many  flow  problems  [1].  The 
reduction  of  wind-timnel  occupancy  is  only  a  direct  benefit  of  CFD.  Detailed 
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flow  information  that  is  unattainable  from  the  wind-tuimel  test  often  can  be 
obtained  by  numerical  simulations.  This  information  about  the  underlying  flows 
can  lead  to  a  better  understanding  and  ultimately  to  a  better  design.  In  fact,  the 
implementation  of  CFD  has  fostered  a  revolution  in  the  design  process  of  flight 
vehicles. 

Current  CFD  methods  have  only  demonstrated  an  ability  to  simulate  flows 
about  complex  geometries  with  simple  physics  or  about  simple  geometries  with 
more  complex  physics.  The  panel  methods  are  unique  in  that  they  provide  a 
capabihty  for  solving  the  flow  about  completely  general  configuration.  Their 
principal  limitation  is  that  they  are  restricted  to  simple  physics  as  modeled  by  the 
linear  equation.  The  Euler  code  provides  more  flow  information  than  the  panel 
method;  however,  the  important  viscous  effects  are  excluded  in  the  approach. 
For  more  realistic  flow  simulations,  the  application  of  Navier-Stokes  equations  is 
required.  However,  the  use  of  full  Navier-Stokes  equations  is  simply  too 
expensive  to  be  practical  for  today's  computers.  An  alternative  is  to  use  the  thin- 
layer  Navier-Stokes  equations.  In  general,  the  governing  equations  are 
approximated  by  a  numerical  scheme  and  then  solved  in  either  a  structured  or  an 
unstructured  grid  network.  The  latter  is  quite  flexible  in  gridding  but  it  requires 
more  work  in  bookkeeping  and  computer  coding.  In  the  structured  grid 
approach,  the  governing  equations  are  first  approximated  by  either  a  finite- 
difference  or  a  finite-volume  scheme,  and  then  solved  in  a  boimdary-fitted 
curvilinear  coordinate  system.  This  approach  has  advantages  in  boundary 
condition  treatment  and  computer  coding;  however,  it  may  have  difficulty  in 
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generating  good  grid  system  for  complex  configurations.  As  evidenced  by  papers 
published  in  journals  and  presented  at  conferences,  e.g.  [2-4],  the  flow  simulations 
involving  simple  geometries  have  progressed  considerably;  however,  the  capability 
for  turbulent  flow  problems  involving  complex  3-D  configurations  is  still  in  a 
development  state.  Hence,  further  research  and  development  effort  is  needed  to 
advance  the  CFD  methods  for  more  complex  geometry  flow  simulations. 

For  flow  past  a  complex  configuration,  such  as  a  flight  vehicle,  the 
generation  of  an  acceptable  single  structured  grid  for  the  entire  flow  field  is 
rather  difficult.  It  is  much  easier  to  divide  the  flow  field  into  several  subregions, 
termed  blocks,  and  then  generate  a  good  grid  for  each  subregion.  Hence,  it  is 
natural  to  develop  a  solution  method  based  on  the  multi-block  grid  network. 
There  have  been  a  number  of  multi-block  solution  techniques  proposed  for  the 
Euler  simulation  of  aircraft  models  and  other  complex  configurations  [5-7]. 
However,  the  topology  of  block  gridding  for  viscous  turbulent  flow  simulation  is 
more  involved  than  that  for  inviscid  flows  and  requires  a  lot  more  grid  points  to 
resolve  the  thin  viscous  layer.  In  general,  a  good  block  grid  topology  for  Euler 
simulation  may  not  be  a  proper  one  for  Navier-Stokes  sunulation. 

There  are  also  zonal  methods  developed  to  use  different  equation  sets  in 
different  regions  of  the  flow  field.  In  general,  the  flow  field  is  divided  into 
inviscid  and  viscous  zones,  and  then  the  Euler  equations  are  appUed  in  the 
inviscid  zone,  while  the  boundary-layer  or  Navier-Stokes  equations  are  used  in  the 
viscous  zone.  A  zonal  approach  that  was  based  on  the  Euler/thin-layer  Navier- 
Stokes  equations  was  proposed  and  tested  successfully  for  transonic  wing  flow  [8]. 


Another  zonal  algorithm  that  applied  the  Euler  and  thin-layer  Navier-Stokes 
equations  for  simulating  the  flow  field  of  isolated  wings  was  reported  in  [9].  The 
flow  field  was  divided  into  four  zones;  two  inviscid  zones  with  coarse  grids  and 
two  viscous  zones  with  clustered  grids.  The  computed  results  were  in  good 
agreement  with  experimental  data  for  the  surface  pressure,  except  in  the 
immediate  vicinity  of  the  tip  and  in  the  shock-induced  separated  region.  The 
zonal  approach  also  has  been  used  to  simulate  flow  problems  for  more  complex 
geometry.  A  transonic  viscous  flow  past  an  F-16A  wing-fuselage  configuration  has 
been  simulated  by  a  zonal  approach  [10].  The  flow  field  was  divided  into  as 
many  as  16  zones;  in  the  inner  zones  adjacent  to  no-slip  surfaces  the  thin-layer 
Navier-Stokes  equations  are  solved,  while  in  the  outer  zones  the  Euler  equations 
are  used.  The  prediction  of  wing  surface  pressures  was  quite  good,  except  at  the 
leading  edge  and  the  aft-shock  position.  The  zonal  approach  has  been  shown  to 
be  rather  efficient.  However,  the  difficulties  are  to  locate  properly  the  interfaces 
for  the  inviscid  and  viscous  zones  and  to  patch  the  solutions  that  are  obtained 
from  different  equation  sets  at  these  interfaces.  Recently,  a  zonal  method  with 
Chimera  overset  grid  scheme  was  investigated  for  subsonic  turbulent  flow  about 
the  F-18  fuselage  forebody  and  the  combined  wing-fuselage  [11].  Special  coding 
was  required  for  the  overlapped  grid  region,  but  the  computed  results  have  been 
shown  to  be  in  good  agreement  with  flight-test  flow  visualization  and  surface 
pressure  measurements. 

In  this  study,  a  novel  multi-block  computational  method  is  proposed  and 
explored  for  the  simulation  of  transonic  turbulent  flows  over  complex 


configurations.  In  this  method,  the  flow  domain  is  divided  into  several  contiguous 
blocks  such  that  each  block  is  partially  boimded  by  a  solid  surface.  For  ease  in 
applying  thin-layer  approximation  and  the  Baldwin-Lomax  turbulence  model,  the 
sohd  smface  is  mapped  onto  an  entire  boundary  plane  in  computational  space. 
In  the  solution  process,  each  block  is  then  treated  as  an  independent  flow 
problem  governed  by  the  same  set  of  time  dependent  thin-layer  Navier-Stokes 
equations.  The  interface  boimdary  conditions  between  blocks  are  updated  at 
every  time  step  of  the  solution  algorithm.  The  solution  algorithm  with  distance 
weighted  interpolation  for  updating  the  interface  boundary  condition  and  the 
computer  coding  are  rather  simple  and  straightforward.  The  multi-block  method 
offers  several  distinct  advantages  over  the  single  block  computation.  In  fact,  the 
use  of  block  grids  makes  the  problem  of  simulating  flow  fields  about  complex 
geometry  more  tenable.  If  a  geometric  feature  is  changed,  only  the  related  block 
requires  to  be  modified  without  completely  redoing  the  basic  grid  topology.  Also, 
the  solution  procedures  are  adapted  to  the  modern  computer  architecmre;  in 
parallel  processing,  each  block  can  be  solved  at  different  processors;  or  in  case  * 
short  of  main  memory,  each  block  can  be  solved  at  a  time  while  the  other  blocks 
remain  on  extended  storage. 

A  transonic  turbulent  flow  over  a  wing-fuselage  configuration  is 
investigated  for  the  development  of  the  proposed  method.  The  flow  field  is 
properly  divided  into  six  contiguous  blocks.  A  coarse  base  grid  was  first 
generated  and  used  for  the  flow  simulation,  and  then  a  refined  grid  was 
generated  from  the  base  grid  by  halving  the  grid  spacing  along  the  wing  in  the 
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chordwise  direction.  A  flow  simulation  package  that  includes  grid  generation, 
flow  solver,  and  plotting  program,  was  developed  and  applied  to  the  six-block 
flow  problem.  The  method  was  first  investigated  for  a  transonic  tiu-bulent  flow  of 
Mach  0.8  past  the  wing-fuselage  configuration  at  zero  angle  of  attack,  then 
simulations  for  different  angles  of  attack,  a =4"  and  a =-3°,  were  conducted  to 
gain  further  insight  and  understanding  on  the  solution  algorithm.  With  the  use  of 
the  refined  grid,  the  calculated  wing  surface  pressure  distribution  agrees  well  with 
measured  data.  Nimierical  experiments  conducted  also  show  that  special 
measures  and  care  are  required  in  generating  good  grid  systems  involving 
excessive  distortion  between  the  physical  and  computational  domains;  however, 
the  results  obtained  indicated  that  the  proposed  method  is  a  very  promising 
approach  for  the  simulation  of  complex  configuration  aerodynamics. 


CHAPTER  n 
GOVERNING  EQUATIONS 

Some  theoretical  background  that  required  for  transonic  turbulent  flow 
simulations  are  described  in  this  chapter.  The  Navier-Stokes  equations  and  ideal 
gas  equation  of  state  are  first  discussed.  Then,  they  are  non-dimensionlized  and 
transformed  to  a  curvilinear  coordinate  system  for  ease  in  numerical  appUcations. 
The  thin-layer  Navier-Stokes  equations  valid  for  high  Reynolds  number  flows  are 
obtained  by  neglecting  the  diffusion  terms  along  the  du-ection  of  the  soUd  surface. 
The  time-averaged  thin-layer  Navier-stokes  equations  are  used  for  turbulent  flow 
simulations  and  the  Baldwin-Lomax  model  [12]  is  used  as  the  turbulence  closure. 

2.1  Navier-Stokes  Equations  and  Equation  of  State 

2.1.1  Navier-Stokes  Equations 

The  fundamental  equations  of  fluid  dynamics  are  based  on  the  following 
universal  laws  of  conservation:  (a)  conservation  of  mass;  (b)  conservation  of 
momentum;  (c)  conservation  of  energy.  The  Navier-Stokes  equations  that  govern 
most  of  the  compressible  flow  problems  can  be  obtained  by  applying  those 
universal  laws  to  a  fluid  flow  [13-14]. 
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Continuity  equation.  The  conservation  of  mass  law  applied  to  a  fluid 
passing  through  an  infinitesimal,  fixed  control  volume  without  sources  or  sinks  of 
mass  yields 

dp      dpyx     dps      apw      ^  ,  , 

—  +  +  +  =  0  (2.1) 
at     ax      ay      az  ^  ' 

where  p  is  the  fluid  density  and  u,  v,  w  are  the  components  of  fluid  velocity  in  x, 

y,  z  direction,  respectively. 

Momentum  equations.  The  conservation  of  momentum  law  appUed  to  the 

control  volume  while  neglecting  body  forces  reads 


i£^4--^(pu^+p)+-^+i^  = 


zx 


at     ax        "     dy       az  ax      ay  az 

at     ax        ay  az  ax      ay  az 

a^w   apwu     apwv     a   ,   ,    ^  8t„     dr  ar 

 +—  +—  +   (pw^  +  p)  =   S.+  .L_S.+  "'^a 


(2.2) 


at     ax       ay       az  ax      ay  az 

where  the  components  of  the  viscous  stress  tensor  ry  for  a  Newtonian  fluid  are 
given  by  the  constitutive  equation 

au     3v     aw.  ^   au     2     au    av  aw 


'"xx 


=  A(  — +— +— )  +  2m—  =  ±^(2^-^.^) 
ax     ay     az        ax     3  ^  ax    ay    az  ^ 


"     'ay     sz     5x  '    '  ay     3  '^^  ay    az    ax  ' 
_  ,  3w    3u     av        aw     2     aw  au    av  , 


au  av 


r 


av  aw 


,  aw     au  , 


IS 

is  a 


Here,  p  is  the  thermodynamic  pressure  and  fi  is  the  dynamic  viscosity.  The  fluid 
assumed  to  be  isotropic  and  satisfying  the  Stokes  condition,  A  =  -2/3  n,  which 
good  approximation  for  flow  problems  without  relaxation  process. 

For  air,  the  viscosity  ^  depends  mainly  upon  the  temperature  and  can  be 
estimated  by  using  an  interpolation  formula  based  on  Sutherland's  theory  of 
viscosity,  i.e., 

,  T,L5  T„  +  S 

M  =  f^.(^r         g  (2.2b) 

where     denotes  the  viscosity  at  the  reference  temperature  T„,  and  S  is  a 
constant  which  assimies  the  value  198.6"R  [15].  '  '  - 

Energy  equation.  The  conservation  of  energy  law  applied  to  the  control 
volume  and  assuming  no  external  heat  addition  results 

-^[(E+p)u]+  i-[(E+p)v]+  4-[(E+p)w]  =  ^  (2.3) 

01     ax  ay  az  ax     ay  az 

where  ■!  . 

aT 

^,  =  Ur„  +  Vr^  +  Wr„+k—  ^  .  ■  - 

aT    ''■  .■ 

^  =  Ury^  +  VTyy  +  Wry,  +k—  "  -     ,  (2.3a) 

aT 

=  Ur„+Vr^  +  Wr„  +k—  , 
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Here,  E  is  the  total  energy  per  unit  volume,  T  is  the  temperature  and  k  is  the 
thermal  conductivity.  The  heat  loss  by  conduction  qj  is  related  to  the 
temperature  gradient  by  Fourier's  law  of  conduction 
aT 

=  -k—  (2.4) 

The  equations  described  above,  Eqs.(2. 1-2.3),  constitute  a  complete  set  of 
the  time  dependent  Navier-Stokes  equations  that  is  valid  for  a  class  of 
compressible  flow  problems  without  heat  addition  and  chemical  reaction  (or 
relaxation).  The  Navier-Stokes  equations  includes  one  continuity  equation,  three 
momentum  equations  and  one  energy  equation,  but  there  are  seven  unknowns, 
namely,  p,  u,  v,  w,  p,  E,  and  T.  The  relationships  between  the  thermodynamic 
variables  (p,  p,  T,  E)  are  established  by  the  equation  of  state  of  the  working  fluid. 

2.1.2  Equation  of  State 

In  this  study,  air  is  the  assumed  working  fluid  and  the  perfect  gas  equation 
of  state  is  employed.  For  a  perfect  gas,  the  thermal  and  caloric  equations  of  state 
are 

p  =  pRT,  e  =  QT    (  Q  =  constant  )  (2.5) 
respectively.  Here,  R  is  the  universal  gas  constant,  e  is  the  internal  energy. 
Some  other  useful  relations  are 

h  =  CpT,  7  =  Cp/Q  a^  =  -yRT  (2.6) 
where  h  is  the  enthalpy,  Cp  and  Q  are  the  specific  heat  at  constant  pressure  and 
volume,  respectively,  7  is  the  ratio  of  specific  heats,  and  a  is  the  local  speed  of 
sound. 


If  only  internal  energy  e  and  kinetic  energy  are  considered  significant,  the 
total  energy  E  can  be  written  as 

E  =  p[e  +  (u^+v^+w^)/2] 
Consequently,  one  obtains  two  specific  equations  relating  p,  T  to  the  dependent 
variables  p,  u,  v,  w,  and  E,  i.e., 

1     c  1  (2.7) 

If  expresses  heat  conduction  qj  in  terms  of  internal  energy  through  the  use  of 
caloric  equation  of  state,  one  can  rewrite  Eq.(2.4)  as  follows. 
,  aT        k    ae        yu  ae 

Here,  the  Prandtl  number  Pr  is  defined  as  Pr  =  C^ti/k,  for  air  Pr  =  0.72, 
approximately,  and  the  internal  energy  is  given  by 

e  =  ^-yCu'+v^+w^)  (2.8a) 
and  the  /9's  in  Eq.(2.3a)  can  also  be  rewritten  as 

iii  de 

=  Ur„  +  Vr„  +  Wr„  + 


"     Pr  ax 

ae 

=  Ur^  +  Vr^+Wr^  +  —  —  (2.9) 


=  Ur„  +  Vr„  +  Wr  + 


Pr  ay 
7/i  ae 


"      °     Pr  az 


2.1,3  Vector  Form 
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The  Navier-Stokes  equations,  Eqs.(2. 1-2.3),  form  a  system  of  five  equations 
and  can  be  cast  in  vector  form  as 


aq     aE     aF     aG      aEL    aF,  aG, 

at     ax     ay     az       ax     ay  az 


(2.10) 


where  q  is  an  unknown  vector,  E,  F,  G  are  the  flux  vectors  and  E^  F^  G^  are  the 
viscous  flux  vectors.  The  explicit  expressions  for  these  vectors  are 


q  = 


p 

\  1 

pU 

pu  +p 

'  pV 

E  =  < 

pUV 

,      F  = 

pW 

puw 

[e  J 

L(E+p)uJ 

v 

pV 
pVU 

pv^+p 

pVW 


G  = 


pW 

pWU 

pWV 

3W^  + 


pW^  +  P 

Ue+p)wJ 


0 


xz 
I9x  J 


F  = 


'0  " 


ro 


(2.10a) 


'  T 


zy 


[fir) 


in  which  ry  are  given  by  Eq.(2.2a),  ^'s  are  defined  in  Eq.(2.3a),  and  p,  T  are 
related  to  q  by  Eq.(2.7). 

2.2  Non-dimensionalization 


To  reduce  the  number  of  parameters  that  characterize  the  same  class  of 
problems,  a  dimensional  analysis  is  employed.  The  non-dimensional  quantities 
are  obtained  by  choosing  characteristic  scales  in  the  following  manner: 


where  the  non-dimensional  variables  are  denoted  by  an  asterisk  *,  the  freestream 
conditions  by  «,  a„  is  the  freestream  speed  of  sound,  and  L  is  a  characteristic 
length  scale.  "     ■      '  - 

By  substituting  the  above  relations  into  the  governing  equations,  Eq.(2.10), 
the  non-dimensional  equations  in  vector  form  are 
aq*  ,  aE*    aF*    aG*    „     dE*  dF*   aG ' 

-T-  +  — s-  +  — s-  +  — r  =  Re  X  — ^  +  — ^  +  — /)  (2.12) 

at      ax     ay     az  ^  ax     ay     az  ^  .  ;  ^  ^ 

where  the  Reynolds  number,  Re,  is  defined  as  v. 
Re  =  ■ 

and  the  relationships  from  the  perfect  gas  equation  of  state,  Eq,(2.7)  become 

p*  =  (^-l)[E*Ap\u*^+v-2+w-2)]=  (-y-l)pV     C  V  • 

E*  1  ■  '  " 

1^  =  7(7-l)[— -^(u'^+v'^+w'^)]  ''^^     ^   •  (2.13) 

The  unknown  vector  q*  and  flux  vectors  E*,  F*,  G*,  E/,  F/  and  G*  have  the  same 
form  as  in  Eq.(2.10a),  except  that  all  the  quantities  are  now  non-dimensional. 

In  summary,  the  dimensionless  equations  govern  a  class  of  flow  problem. 
Only  two  parameters,  namely,  the  Reynolds  number,  Re,  and  the  Prandtl  number, 
Pr,  appear  in  the  equations,  the  Mach  number  is  dropped  out  due  to  the  use  of 
a,  as  characteristic  velocity  scale.  In  the  following,  the  asterisks  *  are  removed 
from  the  non-dimensional  equations  except  where  noted. 
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2.3  Transformed  Navier-Stokes  Equations 


In  the  finite-difference  approach,  to  enhance  numerical  accvu-acy  and 
efficiency  in  applying  boundary  conditions,  the  governing  equations  are 
transformed  from  the  Cartesian  coordinates  to  a  boundary-fitted  curvilinear 
coordinate  system.  Consequently,  a  general  computational  code  can  be  developed 
such  that  the  computation  is  performed  in  an  uniformly-discretized  computational 
domain  [16]. 

If  one  introduces  a  generalized  coordinate  system 
e  =  ?(  x,y,z,t ) 

n  =  r,(  x,y,z,t )  (2.14) 

r  =  r(  x,y,z,t ) 

r  =  t 

and  applies  the  chain  rule  of  partial  differentiation,  the  non-dimensional 
governing  equation,  Eq.(2.12),  becomes 

q.  +  QeC.  +  q.-^. + q^  r .  +  (C^^  +  v^fi,  +  ^^,)  +  (CyF^  +  VyF,  +  TyF,)  +  (e,G^  +  r,fi^  +  r,G,) 
=  Re-^[(^  A,  +  r,^^  +  rxE^)  +        +  VyF^  +  CyF^)  +  (C,G^  +  vfi^  +  r,G^)]    (2. 15) 
The  metric  coefficients  appearing  in  these  equations  can  be  determined  by 
the  matrix  relation  as  follow 


i^  ' 

^?  ^f) 

Vy  »?, 

y?  y.  Vc  yr 

k 

'<}  . 

T.f  Z„  Zf. 

.  o'  0"     \  . 

or 
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f/x 

'y 


X  ''y  ''z 
■  fx  f y  f z  . 


Xj. 

y?  y.  yf 

Zf 


-1 


=  J 


■(yfZf-yfZe)    /fZ^-x^z  -(x^y^-x,y^) 
■(XfZ^-x,Zj)     x^y^-x^y^  J 


yez.-y.Zf 


'^x  Cy 

^0 

' »»» 

'?x  '?y 

'?Z 

"yr ' 

I  Tx  Ty 

fz  J 

Xj  X^  Xj. 

y?  y.  yf 

z. 


where  J  is  Jacobian  of  the  transformation  and  J'\  representing  the  volume  of  the 
grid  cell  in  physical  space,  is  given  by 
_  a(x,y,z) 
a(e,i7.f) 

=  Xe(y.ZryfZ,)-x,(yeZf-yfZ^)+x^(y^z^-y^Zj)  .  .  (2.17) 

In  order  to  put  the  transformed  equations,  Eq.(2,15),  into  conservation-law 
form,  the  equation  is  first  divided  by  the  Jacobian  and  then  rearranged  to  give 
(qJ/J  +  (l.qf  +  ?xEe  +  eyF,  +  e.G,)/J  +  („,q^  +  ^  JE,  +       +  ^^g,)/j 
+  (r.qf  +  rxE,  +  fyF,  +  f  ,G,)/J  =  Re-^[(e      +  e/,,  +  ?,G,,) 

and  then  rewritten  into  the  following  conservation  law  form 

(q/J)r  +  (a.q+e,E+eyF+C,G)/J)^  +  (M+»7^+r,yF+r,,G)/J)„ 

+  ((r.q+fxE+fyF+r,G)/J),  =  Re-^[((eA+^yF,+  ?,GJ/J),  +  . 

+  ((.7A+'7yF,+  r,,G,)/J),  +  ((rJE,+  ryF,+  f,G,)/J),]       ,  (2.18) 

in  which  the  follov^ing  metric  identities  [10]  have  been  apphed. 
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(2.19) 


— (7^)+— (f)  =  0 
a?  ^  J  ^    at;  M      ar  M  ^ 

The  transformed  equation,  Eq,(2.18),  has  the  same  form  as  the  original  equation, 

Eq.(2. 12),  and  can  be  written  as    -    v  ■       V  '  .v 


ar     ae     a»7     ar  a$     a?;  ac 


(2.20) 


where 


q  = 

q/J, 

E 

=  (^,q+|^+C/+|,G)/J, 

F  = 

(,7.q+»7^  +  r7yF+r,,G)/J, 

A 

G 

=  (r,q+fJE+fyF+f,G)/J,  (2.20a) 

A 

=  (^A+?,Ft+?zG.)/J, 

A 

=  („A+'/yFT+'7zG,)/J, 

=  (rA+ryF,+f,G,)/j 

The  explicit  expressions  for  the  transformed  unknown  vector,  flux  vectors 
and  viscous  flux  vectors  are 


-  T-i 


q  =  J 


p 

pV 
pVi 

E 


F  =  J-^ 


E  =  J 


-1 


/juU+^j) 
pvU+^p 
pwU+^j) 
L(E+p)U-^j5j 


puV+r/J) 
pvV+»7yP 
pwV+r/J) 

(E+p)V-„j)J 
0 

Cx'-xy+Vyy  +  ^z''xy 

^x'-xt+Vp'^^^^^ 
^x^x+?A+^A 


G  =  J-^ 


puW+fj) 
pvW+Typ 
pwW+fj) 

[(E+p)w-rj)J 


F,  =  J-^ 


0 


'?x'"xx+Vyx+''z'"zx 
'7x^xy+Vyy  +  ''z'-zy 


(2.20b) 


0 

where  U,  V  and  W  are  contravariant  velocity  components  given  by 

V  =  f7,+T7,u+r;yV+»7,w  (2.20c) 

w  =  r,+f,u+fyV+r^w  ' 

It  is  noted  that  the  velocity  gradients  in  Ty,  Eq.  (2.2a),  and  the 
temperature  gradient  in  /9's,  Eq.(2.3a),  have  to  be  transformed  into  computational 
space  also,  e.g.,   •  -k 

ax     ay     az  '  ax 


-  ^r^  3u       au       au         av^     av       av  ^ 


a?  a»; 


ar 


.  aw  ^  _  aw       aw       ^     au       au  au 
+  (— — »7z+  — rz)]+24(— -Ffx)] 


ae 


a»7 


ar 


a^ 


at? 


ar 


2.4  Thin-Layer  Approximation 


The  numerical  simulation  of  complete  Navier-Stokes  equations,  Eq.(2.20), 
is  in  general  rather  time  consuming  and  needs  a  huge  amount  of  computer 
storage.  For  high  Reynolds  number  flows,  however,  the  viscous  effects  are  often 
confined  to  a  thin  layer  near  no-slip  boundaries,  called  the  boundary-layer, 
outside  of  which  the  flow  can  be  considered  to  be  mostly  inviscid.  With  this  in 
mind,  it  is  reasonable  to  assume  that  the  diffusion  processes  along  the  body 
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surface  are  negligibly  small  in  comparison  with  those  in  the  direction  normal  to 
the  body  surface.  As  a  consequence  of  this  assumption,  all  viscous  terms 
containing  derivatives  parallel  to  the  body  surface  are  dropped  since  they  are 
substantially  smaller  than  viscous  terms  containing  derivatives  normal  to  the 
surface,  while  all  other  terms  in  the  momentum  equations  are  retained.  One  of 
the  advantages  of  retaining  the  terms  which  are  normally  neglected  in  boundary- 
layer  theory  is  that  separated  and  reverse  flow  regions  can  also  be  computed. 
Also,  flows  which  contain  a  large  normal  pressure  gradient  can  be  readily 
computed. 

Another  justification  of  the  use  of  a  thin-layer  approximation  can  be  given 
from  a  computational  viewpoint.  In  practice,  a  substantial  fraction  of  the 
available  computer  storage  and  time  is  expended  in  resolving  the  cross-stream 
gradients  in  the  boundary  layer  since  a  highly  stretched  grid  is  required,  and  the 
resolution  along  the  body  is  treated  as  what  is  needed  in  inviscid  flow.  As  a 
result,  even  though  the  full  derivatives  are  retained  in  the  equations,  the  gradients 
along  the  body  are  not  resolved  in  an  adequate  manner.  Hence,  one  can  drop 
those  terms  that  are  not  being  adequately  resolved,  provided  that  they  are 
reasonably  small. 

Since  the  thin-layer  model  requires  a  boundary-layer  type  of  coordinate 
system,  it  is  assumed  that  the  sohd  body  surface  is  mapped  onto  an  entire  C»7- 
plane  in  the  transformed  space,  then  the  thin-layer  approximation  requires  that 
the  e  and  f;  derivatives  in  the  viscous  terms  be  dropped.  Upon  simplifying  the 
complete  Navier-Stokes  equations,  Eq.(2.20),  the  system  of  governing  equations 
becomes 

i 


aq  ,  aE  ,  aF  ,  aG  _  aS 


dr 

where 


ae  dri 


=  Re-\— ) 
ac 


(2.21) 


S  =  J-^ 


0  * 

miUj.  +  m25-, 
miVf  +  mzfy 

I  mim3  +  m2(r^u  +  ryV+r,w)  J 


(2.21a) 


^2  =  M(r,Uf  +  fyVf  +  r,w^)/3  (2.21b) 

mj  =  (u2+v2+w2)^/2  +  PT\y-ir\r), 

Note  that  the  notation    has  been  removed  from  the  above  equations  for 

simplicity  and  the  heat  conduction  term  are  replaced  by 

jfi    ae    _       /i  aT 
'^""Pr   ar    ~  "(7-l)Pr  "aT  ^^'^^^ 

The  thin-layer  Navier-Stokes  equations,  Eq.(2.21),  in  conjunction  with  the 

equation  of  state,  Eq.(2.13),  constitute  the  required  governing  equations  for  this 

study.  They  are  valid  for  both  laminar  and  turbulent  flow  problems.  However, 

due  to  the  widely  disparate  scales  of  turbulence,  it  is  impractical  to  apply  the 

equations  directly  for  turbulent  flow  simulation.  An  alternative  is  to  use  the 

time-averaged  equations  coupled  with  a  turbulence  model. 


2.5  Turbulence  Model 


The  time-averaged  equation  is  obtained  by  decomposing  the  dependent 
variables  into  time  mean  and  fluctuating  components  and  then  time-averaging  the 
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entire  equation.  In  the  averaging  process,  new  unknowns  will  be  introduced  and 
must  be  treated  by  some  approximating  techniques  so  as  to  close  the  problem.  In 
this  study,  the  mass-weighted  time-averaged  equations  are  used  for  turbulence 
flow  simulation  with  the  eddy  viscosity  that  appears  in  the  averaged  equations 
provided  by  the  Baldwin-Lomax  two-layer  algebraic  model.  In  this  approch,  the 
mass-weighted  time-averaged  thin-layer  Navier-Stokes  equations  will  have  the 
same  form  as  the  original  equations,  Eq.(2.21),  while  all  the  dependent  variables 
are  time-averaged.  The  effective  coefficient  of  viscosity  n  is  comprised  of  two 
parts,  the  laminar  viscosity  /i,  and  the  tiu-bulent  eddy  viscosity  /i,.  The  laminar 
viscosity  A*,  is  obtained  from  Sutherland's  theory  of  viscosity,  Eq.(2.2b),  and  the 
eddy  viscosity  ^l^  is  provided  by  a  turbulence  model  [16]. 

Similarly,  the  effective  thermal  conductivity  k  is  also  comprised  of  two 
parts,  k,  and  k,.  The  heat  conduction  is  expressed  in  dimensional  quantities  as 

The  heat  conduction  term  in  Eq.(2.22)  is  then  replaced  by  its  dimensionless 
counterpart 

Here,  the  turbulent  Prandtl  number  Pr,  is  defined  as  Pr,  =  ^x^c^/k^,  and  is  about 
0.9  for  air  [10]. 

In  this  approach,  the  eddy  viscosity  /i,  is  the  only  quantity  that  has  to  be 
modeled.  In  what  follows,  the  theoretical  background  about  the  algebraic  model 
will  be  discussed  first,  and  then  the  Baldwin  Lomax  model  is  described. 
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2.5.1  Two  layer  algebraic  model 

In  analogy  with  the  coefficient  of  viscosity  in  Stoke's  law  for  laminar  flow, 
Boussinesq  [17]  introduced  a  mixing  coefficient,  /i„  or  referred  to  as  eddy- 
viscosity,  for  the  Reynolds-stress  that  appears  in  the  time-averaged  equation, 
au 

^  =  /^.-T—  (2.24) 
ay 

where  u  is  the  streamwise  velocity,  y  is  the  spatial  coordinate  normal  to  the  solid 
siuface,  and  3u/ay  is  assumed  to  be  the  only  significant  component  of  strain  rate. 

With  Prandtl's  mixing  length  hypothesis,  the  coefficient  n^,  is  then  related  to 
mean  velocity  field  by  'v- 


du 
dy 


.    .         .   :r  ■■  (2.25) 

where  1  is  the  mixing  length.  ' 

The  algebraic  model  is  then  to  establish  a  relation  between  the  mixing  length 
and  the  characteristic  length  of  the  respective  flow.  For  simple  wall  shear  flows, 
the  mixing  length  in  the  law  of  the  wall  region  is  proportional  to  the  normal 
distance  from  the  wall,  i.e.,  1  =  ky,  where  k  =  0.41  is  a  universal  constant 
determined  experimentally.  In  the  viscous  sublayer,  however,  the  mixing  length 
theory  is  not  valid.  A  corrected  form  of  the  mixing  length  valid  down  to  the  wall 
was  deduced  by  van  Driest.  He  introduced  a  damping  factor  and  wrote 

1  =  ky[l-exp(-yVA^)]  .     ,  (2.26) 
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where  y"*"  =  ypujn  is  the  law  of  the  wall  coordinate,  A"^  denotes  an  empirical 
constant  which  equals  26  for  flat-plate  flow,  and  u^,  the  friction  velocity,  is 
defined  as      =  {rjpf',  where  r„  is  the  shear  stress  on  the  wall. 

In  the  law  of  the  wake  region,  the  Clauser  model  is  commonly  used  in 
conjunction  with  a  factor  to  account  for  the  intermittent  effect.  The  eddy 
viscosity  is  then  written  as 


where  the  value     is  obtamed  experimentally  and  assumed  to  be  0.0168,  U  is  the 
streamvdse  velocity  on  the  edge  of  the  boundary  layer,  5*  is  the  boundary  layer 
thickness,  7  is  the  intermittency  factor  deduced  from  a  curve  fit  of  measured  data 
and  given  by  Reynolds  and  Cebeci  [18]  as 


where  s  is  the  boimdary  layer  thickness. 

Based  on  the  above  arguments,  Cebeci  and  Smith  [19]  proposed  a  well 
known  two-layer  algebraic  model  as  follows. 
In  the  iimer-layer: 


(2.27) 


7  =  [l  +  5.5(y/5)'] 


(2.27a) 


(2.28) 


1  =  ky[l-exp(-yVA*)] 


(228a) 


In  the  outer-layer: 


A*,  =  k,U5*7 


(2.28b) 


7  =  [l+5.5(yA)«] 


(2.28c) 


This  model  has  been  shown  to  be  in  excellent  agreement  with  experimental    ■  ' ; , 
observations  for  simple  shear  flows.  However,  the  model  requires  the  boundary 
layer  thickness  and  the  edge  velocity  which  constitute  a  practical  disadvantage  in 
the  implementation  of  the  model.  In  an  attempt  to  overcome  this  difficulty, 
Baldwin  and  Lomax  [12]  modified  the  outer  layer  formulation  to  eliminate  the 
necessity  of  finding  the  boundary  layer  thickness.  In  addition,  they  replaced  the 
velocity  gradient  by  the  vorticity,  which  is  invariant  under  coordinate 
transformation,  so  that  the  model  can  be  easily  applied  to  3-D  problems.  The 
Baldwin-Lomax  model  is  described  next, 

2.5.2  Baldwin-Lomax  Turbulence  model 

The  Baldwin-Lomax  model  is  similar  to  the  Cebeci-Smith  two-layer 
algebraic  turbulence  model.  The  eddy  viscosity  in  a  profile  is  calculated  for  each 
inner  and  outer  layer,  and  then  matching  both  layer  as  follows  [12]: 


where  y,  is  the  smallest  value  of  y  at  which  values  the  inner  and  outer  formulas 
are  equal. 

The  inner  layer  follows  the  Prandtl-Van  Driest  formulation 


where  y  is  the  normal  distance  to  the  surface  and  |w|  is  the  magnimde  of  the 
vorticity 


(2.29) 


(/^.)in»er  =  pl'kl     with  1  =  ky[l-exp(-y  VA*)] 


(2.30) 


(2.31) 
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=  y(py,^J^/f^  is  the  law  of  the  wall  coordinate,  A"^  =  26,  and  /x^  are 

the  local  density,  shear  stress,  and  laminar  viscosity  evaluated  at  wall,  respectively. 
The  eddy  viscosity  for  the  outer  layer  is  given  by 

(M,)outer  =  pk,QpIVakeFkleb(y)       '         -         '  '  •       '  "  (2.32) 

where     =  0.0168  is  the  Clauser  constant  and      =  1.6.  *  .  H  '  , 

The  parameter  F^^^^  is  given  by 

Fwak.  =  mini  J^rf        /p      where       =  0.25     -  -  -  (232a) 

and 

Udif  =  (u'+v^+w^)^^  -  (u^+v^+w^)^^  .  (2.32b) 

The  quantities  y^  and        are  determined  from  the  function 

F(y)  =  yk|[l-exp(-yVA*)]  (2.32c) 
where  is  the  maximum  value  of  F(y)  and  y^  is  the  value  of  y  at  which  F^ 
occurs.  The  Klebanoff  intermittency  factor  Fu,,,(y)  is  given  by 

Fkid.(y)  =  [l  +  5.5(Cu^/y^']-»      and  C^.^  =  0.3  (2.32d) 
The  Baldwin-Lomax  model  has  been  investigated  and  apphed  to  a  variety 
of  2-D  and  3-D  flow  calculations  with  satisfactory  results  [10].  The  model  is 
capable  of  handling  the  compressible  flow  problems  with  modest  flow  separation 
[20].  The  implementation  of  algebraic  models  requires  only  a  small  amount  of 
computer  time  and  storage;  this  is  particularly  important  in  3-D  computations. 


CHAPTER  in 
NUMERICAL  ALGORITHM 


In  numerical  simulations,  the  transformed  governing  equations  described  in 
the  previous  chapter  are  approximated  by  finite-difference  equations  and  then 
solved  in  a  equally  spaced  computational  space.  The  basic  computer  code  used 
in  this  study  is  a  TVD  thin-layer  Navier-Stokes  code  which  is  based  on  an  original 
version  of  thin-layer  Navier-Stokes  code,  ARC3D  of  NASA  Ames  Research 
Center.    The  original  ARC3D  was  based  on  the  Beam  and  Warming  scheme 
[21],  and  has  been  improved  for  its  efficiency,  accuracy  and  convergence  [22];  for 
instance,  there  are  options  for  using  different  numerical  dissipations  for  better 
stability  criteria.  Alternatively,  to  improve  the  robustness  of  the  code,  a 
symmetric  Total  Variation  Diminishmg  (TVD)  scheme  was  implemented  into  the 
original  ARC3D.  This  modified  ARC3D  has  shown  to  be  more  efficient  and 
robust  than  the  original  one  [23-26],  and  hence  is  adopted  in  this  study.  _ 

The  Beam- Warming  scheme  and  the  TVD  scheme  were  constructed  by 
totally  different  philosophies;  however,  in  applications  they  can  be  viewed  as  the 
same  scheme  with  different  numerical  dissipation  terms  [27].  Both  schemes  will 
be  described  and  discussed  in  the  following.  ' 
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3.1  Beam  and  Warming  Scheme 

The  Beam- Wanning  scheme  is  an  non-iterative,  implicit  approximate 
factorization  finite-difference  scheme.  The  implicit  operator  is  factorized 
(spatially  split)  to  obtain  an  ADI-like  (  Alternating  Direction  ImpUcit )  algorithm 
such  that  the  three-dimensional  inversion  process  is  reduced  to  three  one- 
dimensional  processes.  The  basic  algorithm  is  second-order  accurate  in  both  time 
and  space  [28].  The  implementation  of  this  scheme  to  the  system  of  governing 
equations,  Eq.(2.21),  is  described  as  follows. 

3.1.1  Temporal-Differencing 

To  approximate  the  time  derivative  of  Eq.(2.21),  the  well-known  second- 
order  Trapezoidal  Rule  or  Crank-Nicolson  Formula  is  employed  to  give 

^  =  +  3-7^)  +  0[(Ar)2]     where     Aq-=q-*V  (3.1) 

If  one  rewrites  Eq.(2.21)  as 

and  then  substitutes  it  into  Eq.(3.1),  the  following  equation  results 
2  "  3{       3,       ac      ^  3f 

In  the  above  equation,  the  flux  vectors  E,  F,  G,  and  viscous  flux  vector  S  at 
the  advanced  n+ 1  time  step  are  in  general  nonlinear  functions  of  the  unknown 


;  if 
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4 

vector  q.  They  are  linearlized  while  maintaining  second-order  time  accuracy  as 
follows. 

The  flux  vectors  are  linearlized  by  using  a  Taylor  series  expansion  about  q" 
as  follows 

E"+^  =  E"+A"Aq"  +  0[(Aq)^]  where 
pn+i  ^  pn^gn^qn  ^  o[(Aq)^]  wherc 

G°**  =  G"+C"Aq"  +  0[(Aq)^]  where 

The  above  approximation  is  second-order  in  time  since 

Aq"  =  q"-'»-q"  =  0(Ar)    and  0[(Aqf]  =  0[(Arf] 
The  viscous  flux  vector  S  is  linearized  by  using  the  procedure  suggested  by  Steger 
[29].  The  coefficient  of  viscosity  ^  and  thermal  conductivity  k  are  assumed  to  be 
locally  independent  of  Jq  (  recall  that  Jq  is  the  unknown  vector  in  physical 
domain  ),  so  that  elements  of  S  can  be  expressed  in  the  general  form 

Sj  =  J'^aj  -— i       where  a,  is  independent  of  Jq. 

Each  element  is  then  linearized  in  the  following  manner  by  Taylor  series 
expansion 

S,-^  =  S,Vj-l{a,"-^[E(-^r]}jAq/  +  0[(Ar)^] 

Written  in  vector  form,  the  viscous  flux  vector  becomes 

gn+i  ^  s"+J-^M"JAq"  +  0[(Ar)2]  (3,4) 


A  = 


B  = 


C  = 


3E 
aq 

dF 
aq 

aG 
aq 


(3.3) 
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The  matrices  A,  B,  C  and  M  are  so  called  Jacobian  matrices  and  the  elements  of 
these  matrices  are  given  in  Appendix  A. 

After  substituting  all  the  linearized  terms,  Eqs.(3.3-  3.4)  into  Eq.(3.2),  one 

has 

2  ^  3^     dr,     af        ar  ^         ^  ,  ^  .  r ...  ^ 

=  .^,(iE^_aF^_aG_j^^.,^  ,  ^  -  v 

a^     ar;     ar        af  ^  ^  ^  ^ 

where  I  is  the  identity  matrix,  and  Aq"  is  the  delta  form  unknown  to  be  solved 

for.  Note  also  that  a  first-order  time-accurate  Euler  implicit  scheme  can  be 

obtamed  by  simply  replacing  Ar/2  on  the  left-hand  side  of  Eq.(3.5)  by  Ar  [28]. 

3.1.2  Approximate  Factorization 

Solving  Eq.(3.5)  directly  is  very  inefficient  since  it  involves  the  inversion  of 
a  very  large  system  resulted  from  three  dimensions.  Beam  and  Warming  appUed 
an  approximate  factorization  procedure  to  reduce  the  inversion  to  three  one- 
dimensional  inversions.  They  rewrote  Eq.(3.5)  in  the  factored  form 

Since  the  error  introduced  by  this  factorization  procedure  has  the  leading  third- 
order  term 

Ar^._aA_  _aB_^  aB_  _aC     aC  aA  ag" 

4    ac  dr,     dr,  ar     ar  ac 

the  second-order  accuracy  of  the  difference  approximation  is  not  affected. 


3.1.3  Spatial  differencing 

The  spatial  derivatives  appearing  on  the  left-hand  side  of  Eq.(3.6)  are 
approximated  by  second-order  finite-difference  formula,  e.g. 

hence  the  resulting  scheme  possesses  a  block-tridiagonal  structure  that  can  be 
solved  efficiently.  On  the  right-hand  side,  the  same  second-order  differencing  is 
employed  for  the  viscous  term,  while  a  fourth-order  finite  differencing 
approximation,  e.g. 

is  used  for  the  convection  terms  to  keep  the  convection  truncation  error  from 
exceeding  the  magnitude  of  that  of  the  viscous  term,  ^  ;  ■ 

After  applying  spatial  differencing  into  Eq.(3.6),  the  spatial  derivatives  are 
replaced  by  difference  operators,  and  the  finite-difference  form  of  Eq.(2.21)  is 
then  .  ..  "  , 

[1+  y-6,A]"[I+  y^5,B]"[I+  ^(5^C-Re-^S,(J-^MJ))]"Aq- 

=  -ArC^^E+sZ+SfG-Re-^S)"  (3.7) 
This  is  the  basic  Beam-Warming  algorithm.  The  implicit  operator  has 
been  decoupled  to  produce  an  ADI-like  scheme,  and  due  to  the  second-order 
central-difference  operators  employed,  the  algorithm  produces  block  tridiagonal 
systems  for  each  spatial  direction. 
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3.1.4  Numerical  Dissipation 

The  Beam- Wanning  scheme  requires  the  use  of  numerical  dissipations  to 
damp  out  the  spurious  oscillations  that  occurs  when  dealing  with  discontinuities. 
Several  different  numerical  dissipations  have  been  investigated  and  discussed  in 
[30].  Those  employed  in  the  original  ARC3D  code  were  proposed  by  Beam  and 
Warming  [28],  and  Steger  [29]  and  are  described  as  follows: 

The  second-order  implicit  dissipation  operators 

=  -«iJ"\a^J  (3.8) 

are  added  to  the  imphcit  operators  on  the  left-hand  side  of  Eq.(3.7)  in  the  |,  rj 
and  f  directions,  respectively,  and  on  the  right-hand  side  the  following  explicit 
fourth-order  dissipation  term  is  added: 

De  =  -^EJ-'[(VA)'+(V/+(VfA,)']JAq-  (3.9) 
in  which  symbols  a  and  v  are  the  forward  and  backward  difference  operators,  e.g. 

^  =  +  =  -J^i^rki)  +  0[A|] 

and  €j,  eg  are  constant  coefficients  used  to  control  the  amount  of  numerical 
dissipations.  It  was  suggested  to  set  €E=At  and  ei=2€E  so  that  as  the  time  step 
increases,  the  numerical  dissipations  added  relative  to  the  spatial  derivatives  of 
convection  and  diffusion  remains  constant. 


After  inserting  numerical  dissipations  into  Eq.(3.7),  the  finite-difference 
representation  of  Eq.(2.21)  becomes 

[1+  |l-5,A+D,r[I+  ^6^B  +  DJ"[I+  y^(5,C-Re-^5,(J-^MJ))  +  D,]"Aq" 

=  -Ar(SjE+6^F+5j.G  -  Re"^5j.S)"+DE  (3.10) 

3.1.5  Solution  procedure 

For  simplicity,  the  following  notation  is  introduced 


h  =  [I+y-5,A+D,r 


L„  =  [1+  ^«„B  +  D  1" 


Lf  =  [1+  y^(5fC-Re-^5f(J-^MJ))  +  Df]"Aq" 

R"  =  -Ar(5^E+5,,F+5^G-Re-^5^S)"+DE 
where  the  L's  are  operators,  then  the  Beam-Warming  scheme,  Eq.(3.10),  can  be 
written  as 

L,L,L,Aq-  =  R-  (3,11) 
If  one  further  defines  intermediate  solutions  as 
Aq*  =  Lj-Aq" 
Aq    =  L,Aq 

the  unknown  vector  q  at  the  n+1  tune  level  can  be  computed  as  follows. 
LfAq-  =  R» 

L.  • 
/q    =  Aq 

Vq"  =  Aq*  (3  12) 


q"+i  =  q-+Aq" 

The  Beam-Warming  scheme  has  been  shown  to  be  accurate  and  efficient 
for  simple  flow  problems  [10],  however,  for  more  complex  flow  cases  it  can  be 
very  sensitive  to  grid  resolutions  [25]. 

.  ■     :        3.2  TVD  Scheme 

It  is  known  that  finite-difference  schemes  may  produce  oscillations  when 
appUed  across  a  discontinuity.  For  instance,  the  second-order  Lax-Wendroff 
scheme 

u;*i  =     -K  uj^i-u?  )-V2uil-u){  uj"^i-2uj"+u/!i ),  u  =  aAt/Ax  (3.13) 

for  a  linear  hyperbolic  equation 

au  - 
+  a        =0,  a  >  0. 


at  ax 

can  produce  spurious  oscillations  near  the  discontinuity.  The  Total  Variation 
Diminishing  (TVD)  scheme  is  one  of  the  methods  designed  to  eliminate  those 
undesirable  oscillations.  The  notion  of  TVD  was  first  introduced  by  Harten  [31]. 
He  derived  a  set  of  sufficient  conditions  from  which  the  constructed  schemes  will 
not  generate  oscillations  across  discontinuities.  Since  the  theoretical  background 
of  TVD  schemes  is  rather  complex  and  involved,  only  relevant  basic  definitions 
are  given  here.  The  construction  of  a  symmetric  TVD  scheme  and  its  application 
to  thin-layer  Navier-Stokes  equations  is  discussed  in  the  subsequent  sections. 


A  numerical  scheme,  either  explicit  or  implicit,  can  be  written  in  the 
following  operator  form  [32] 

L  u°**  =  R  u"  (3.14) 
where  u""^*  is  the  mesh  function  at  the  advanced  time  step  to  be  solved,  u"  is  the 
known  mesh  function  at  current  time  step,  and  L  and  R  are  some  finite- 
difference  operators  designed  for  the  particular  scheme,  e.g.,  for  the  second-order 
Lax-Wendroff  scheme  given  in  Eq.(3.13) 

L  =  1  ■ 

R  =  1 -.^Vm-MMWh  '  V:l 

Here  and  in  the  following,  Aj+j^  stands  for  the  central  difference  operator,  i.e., 
Aj^.j4U  =  Uj+i  -  Uj.  The  total  variation  of  a  mesh  function  u°  is  defined  to  be 

=iJL  I^Vi-u"!  =jL  ''^ >  '  '  (3.15) 

Here,  the  notation  «  is  used  to  emphasize  that  the  summation  is  taken  over  the 
entire  domain.  If  a  numerical  scheme  applies  to  a  mesh  function  u,  the  total 
variation  of  the  mesh  function  is  non-increasing  as  a  function  of  time,  i.e. 

TV(u"*M<TV(u°)        for  all  n 
the  scheme  is  said  to  be  TVD.  In  practice,  the  following  sufficient  conditions 
derived  by  Harten  [31]  are  served  for  constructing  TVD  schemes: 

TV  (  R  u"   )  <  TV  (  u"  ) 

TV(Lu-^)>TV(u»-^M  (3.16) 
That  is,  if  the  numerical  scheme  shown  in  Eq.(3.13)  satisfies  the  sufficient 
conditions,  Eq.(3.16),  the  scheme  is  TVD. 
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It  has  been  proved  that  a  TVD  scheme  is  monotonicity  preserving  [32], 
that  is,  for  any  monotone  mesh  function  u  that  increasing  (or  decreasing) 
monotonously  with  time,  w  =  Qu  is  also  monotone,  where  Q  is  a  finite-difference 
operator.  Since  a  monotonicity  preserving  scheme  will  not  generate  spurious 
oscillations,  and  neither  will  the  TVD  scheme. 

In  what  follows,  the  construction  of  the  symmetric  TVD  scheme  for  a  one- 
dimension  hyperbolic  equation  will  be  described  first,  then  the  scheme  is  applied 
to  the  Euler  equations,  and  finally  it  is  extended  to  the  thin-layer  Navier-Stokes 
equations. 

3.2.1  Symmetric  TVD  Scheme 

The  symmetric  TVD  scheme  was  developed  by  Yee  [27]  and  it  is  a 
generalization  of  Roe  [33],  Davis  [34]  and  Sweb/s  [35]  TVD  Lax-Wendroff 
scheme.  To  illustrate  the  construction  of  this  scheme,  a  one-dhnension  linear 
hyperbolic  equation  of  the  form 


8u  du 

  +  a  — 

at  ax 


+  a—  =  0  (3.17) 


is  considered.  Sweby  has  shown  that  the  second-order  Lax-Wendroff  scheme  for 
this  equation  can  be  constructed  by  adding  an  anti-diffusive  flux  to  the  first-order 
upwind  scheme  [35],  i.e., 

Uj-**  =       -  i.A..j,u"  -  y2i.(l-..)Aj.j4Aj^^,U»  ,         if  a  >  0. 

uj*»  =     -  ^Aj^^u"  -  Vi«.(l  +  ..)A.^^Aj.j,u-  ,     if  a  <  0.  (3.18) 


where  v  =  aAt/Ax  is  the  Courant  number.  Then  he  introduced  a  limiting 
function  p  and  rewrote  the  above  scheme  as 

i^*i=uJ-.{l  +  y2(l-.)[p(r/)/r/-p(rV,)]}Aj.,^^  ,  if  a  >  0  ^^^^^ 
1^*1=  uJ-u{l-y2(l  +  ..)[p(rj;i)-p(r-p/rj]}Aj,^u-    ,     if  a  <  0 

where  p  is  chosen  to  be  a  function  of  the  ratio  of  consecutive  cell  gradients  rj;  pj 

=  P(rj).      =  ^-mu/Aj+mU,  rj'  =  Aj+^^u/Aj.^^u,  with  the  restrictions  0<p/r<2, 

0<p<2  to  satisfy  the  TVD  constraints,  Eq.(3.16). 

The  Sweb/s  scheme,  Eq.(3.19),  was  reformulated  by  Davis  [34]  as  a  Lax- 

Wendroff  scheme  plus  an  upwind  weighted  numerical  dissipation  term,  i.e. 

[K*(r/)+K-(rj^i)]Aj,^u"-[K*(rV)+K-(r;)]Aj.^u''  (3.20) 


where 


0 
0 


K*(r.-^)=  /  •'(l-'')tl-P(^j')]       if  a  > 
I  0  if  a  <  ( 

KYr-)  ^^-^  (3.20a) 

\'^(1  +  u)\p(t:  )-1]     if  a  <  0 

To  avoid  determining  the  upstream  direction,  Davis  then  modified  the  dissipation 
term  as  follows, 

=  M(i-IH)[i-P(V)]/2 

K-(ij-)  =  IH(l-|H)[l-p(ri-)]/2  .  (3.20b) 

Shortly  after  that.  Roe  [33]  further  generalized  Davis'  scheme,  Eq.(3.20), 
and  suggested  a  new  limiting  function  Q.  Roe's  scheme  takes  the  form 
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**    ^  V 

-  V2|H(i-kl)(i-Qj.H)w"  +  ^2|H(i-M)(i-Qj.h)Vv4u"  (3-21) 

where  the  limiter  Q  depends  on  three  consecutive  gradients  given  by 

Qj+M  =  Q(Aj.WAj+MU  ,  Aj+ii^u/Aj+i^u)  =  Q(r*j  ,  r-_j^i)  (3.21a) 
and 

0<      Qj,H  <2/(l-|H) 

0  <  Qj+^Aj+i  <  2/ 1 1/1  (3.21b) 
0  <  Q^.Jr^j.     <  2/IH 
to  satisfy  TVD  sufficient  conditions,  Eq.(3.16).  There  are  a  variety  of  limiter  Q 
satisfing  the  above  restrictions,  such  as 
Q(r',r'^)  =  max[0,  min(2r",l),  min(r",2)] 

+  max[0,  min(2r*,l),  min(r+,2)]  +  1  (3.21c) 
Roe's  scheme  was  later  modified  and  extended  to  a  one-parameter  family 
TVD  scheme  by  Yee  [27].  She  first  rewrote  Roe'scheme,  Eq.(3.21),  into  the  form 

uy*V  Ae(hP;^^-h-^)  =  uf  +  A(l-e)(hj",^-h;^)  (3.22) 

hj+H  =  [a(Uj+i  +  up-{Aa'Qj^,4+  |a|(l-Qj^H)}Aj+Mu]/2,     a  =  At/Ax  (3.23) 

and  then  simplified  the  numerical  flux  to 

hj+M  =  [a(Uj+i  +  Uj)-|a|(l-Qj^^)Aj^^,u]/2  (3.24) 

Note  that  the  term  Aa^Qj^.^^Aj^.^u/2  was  taken  out  from  Eq.(3.23),  but  it  makes  no 

difference  as  long  as  the  limiter  function  Q  is  chosen  such  that  the  scheme 

satisfies  the  TVD  sufficient  conditions. 
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For  nonlinear  hyperbolic  conservation  law  of  the  form 
au  3f(u) 

the  numerical  flux,  Eq.(3.24),  was  modified  to 

hj+vi  =  [(fj+i+fjH(Vv4)(l-Qj+w)Aj+^u]/2  (3.26) 
where  ^=f(Uj),  Aj^.^4U=Uj+i-Uj,  a(u)  =  df(u)/du  is  characteristic  speed  and 

la(u,)  .VhU  =  0  (^-^ea) 

The  function  *(aj+j4)  that  is  sometimes  referred  to  as  the  coefficient  of  numerical 
viscosity,  is  defined  as 

:;:  *<^)  =  {!^+z^)/2<  ;|z|  J  (3.26b) 

to  avoid  the  entropy  violating  problem  [31],  where  e  is  a  parameter  that  is 
described  in  Ref.[36].  In  order  for  the  above  scheme,  Eqs.(3.22,  3.26),  to  be 
TVD,  that  is,  satisfying  the  TVD  constraints,  Eq.(3.16),  the  following  conditions 
obtained  by  Yee  [37],  must  be  satisfied 
0  <       Qj+H    <  2 

:^  0  <  Qi^Jx\     <  2/[.(l-e)|aj.J]  -  2 

■    0  <  Qj^^/r^^i  <  2/[Kl-e)|aj^m|]-2  ^.  (3.26c) 

Haj^J  <  2/[3(l-e)]  - 
If  one  selects  the  following  Umiter  that  satisfies  the  above  constraints 

Qj+^  =  minmod(Aj.j4U  ,  a^+^u)  +  minmod(Aj+^u  ,  Aj+j^u)  -  Aj^^^u  (3.27) 
then  fi-om  Eqs.(3.26-3.27)  Yee's  numerical  flux  can  be  expressed  as 
,    ^i^^  =  [fj  +  fj+i  +  (3.28) 

fl^j+ii  =  *(VH)(gj+gj+i)  -  2A,^i,u         ^  (3.28a) 
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where  the  minmod  function  niinmod(x,y)  is  defined  as 


minmod(x,y)  =  -  x  ,    x  <  y 


'  0  ,  if  X  y  have  opposite  sign 


ly,    X  >  y 


and 


g.  =  minmod(Aj+v4U  ,  a,.^u) 


(3.28b) 


The  above  scheme,  Eqs.(3.22,  3.28),  form  a  class  of  TVD  schemes  with  the 
numerical  dissipation  terms  centered  and  hence  are  often  referred  to  as 
symmetric  TVD  schemes.  It  is  noted  here  that  the  scheme  is  TVD  for  one- 
dimensional  nonlinear  sccder  hyperbolic  conservation  laws,  and  is  also  TVD  for 
constant  coefficient  hyperbolic  systems.  However,  in  apphcations,  it  is  formally 
extended  to  more  complex  equations,  such  as  the  Euler  equations,  the  thin-layer 
Navier-Stokes  equations,  and  evaluated  by  numerical  experiments. 

3.2.3  Extension  to  Euler  equations 

Extension  of  the  scaler  TVD  scheme  to  systems  of  conservation  laws  is  not 
unique.  Here,  one  follows  the  method  proposed  by  Yee  [38,39].  The  method  is 
accomplished  by  defining  at  each  point  a  "local"  system  of  characteristic  fields  and 
then  applying  the  scheme  to  each  of  the  scalar  characteristic  equations  and  hence 
is  referred  to  as  the  "local  characteristic  approach". 

Consider  the  inviscid  part  of  the  thin-layer  Navier-Stokes  equations, 


Eq.(2.21), 


dq      aE     aF  3G 


=  0 


at     a^     art  a^ 


(3.29) 


Recall  that  the  Jacobian  matrices  of  E,  F,  G  are  A,  B,  and  C,  respectively. 
Denote  the  eigenvalues  of  A,  B,  and  C  as 

(af.a-fja'j.aj.a^),  (a^,a^,a?^,a'^,a^),  (apaji^apa^a^), 
and  Rp  R^,  Rj.  as  the  eigen-matrices  whose  columns  are  eigenvectors  of  A,  B,  C, 
respectively,  and  R"*^  R"*,,  R**j.  as  the  inverses  of  the  corresponding  eigen- 
matrices.  The  explicit  form  of  these  eigenvalues  and  matrices  are  given  in 
Appendix  B. 

In  the  ^-direction,  let  Qj+h^  be  some  symmetric  average  of       and  Qj+i^, 
e.g.  ■     ,  - 

or  Roe's  average 

qj+v4w  =  (^^jW*Jj,k4+^Viw*Jj+iw)/(^\ia+^Viw)  (3-30) 
and  denote         R^^.^^,  R'\+^  as  the  quantities  a'^,  R^,  R"*^  evaluated  at  qj+M,k4- 
The  difference  of  the  characteristic  variables  in  the  local  f -direction  are  defined 
as 

=  RVv4(JjMW<lj+iW-Jik4W/[(Jj+W+W/2]  (3.31) 
One  can  also  define  the  corresponding  parameters  for  the  r)  and  f  directions  in  a 
similar  manner.  ^ 

Using  these  notation  and  applying  Eq.(3.22)  to  each  equations,  a  one- 
parameter  family  of  TVD  scheme  for  the  system  of  equations  can  be  written  as 

qJlM  -  ^(l-e)[(Ej^  wEjW  +  (I^VH4-FjW  +  (<^W.^-^^^^^  (3.32) 
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where  A  =  At  since  A^=Ar7=Ar  =  l 

and  the  numerical  flux  functions  obtained  from  Eq.(3.28)  are  given  by 

^  [Pj,k4''"^j,li+14''"^k+V4*k+M]/2  (3.33) 

the  expression  of  the  elements  of  ^j+j^  are  obtained  from  Eqs.(3.28a-3.28b)  and 
are  given  as 

^j+H  =  *(a'j,J(g'j+g'jM-2a'j^^)  (3.33a) 
gfj     =  minmod(aV.v4,a'j+^)  (3.33b) 
Likewise,  one  can  define  $^+^4, 

3.3.6  Linearization  and  API  decomposition 

In  order  to  solve  Eq.(3.32)  efficiently,  the  left-hand  side  of  Eq.(3.32)  is  first 
linearized  by  the  procedure  described  in  Section  3.1.1.  The  resulting  scheme 
written  in  delta  form  is 
[I  +  Ae[(Hj,^^-Hj.^^  +  (Hj^,H4-Hj^.^^)  +  (H^w^H-Hjw.v,)]Aq» 

= -M(Ej.H,k4-%Mw)  +  (Fj^+M4-FjW  +  (Gjw+H-G^^^^  .  (3.34) 
where  I  is  a  5  X  5  identity  matrix  and  the  H  operators  are  given  by: 

^+MW  ~  tAi+W'""j+K,k4]/2 

Hj.k+M4  =  [Bjj.+M+"j4c+v4.i]/2  (334a) 
in  which 

"j+M,M  =  {R  diag[^'-2*(a')]  R-%^Ai^^ 
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■i    nj^+H4  =      diag[;9'-2*(a')]  R-\^^a^,^  (3.34b) 

"jW+H  =  {R  diag[y3'-2*(a')]  R-'},+mAi+m 
and 

=  *(a^+H)(g^+g^+i)/«j+M      ^    ".    .  J^;^   -  ^  " 

)8k+v4  =  *(ak+v4)(g'k+g'k+i)/«k+v4  "   '  (3-34C) 

Here,  the  expression  diag(z')  denotes  a  diagonal  matrix  with  diagonal  elements  z*. 

To  calculate  the  n's  at  every  time  step  is  quite  costly.  For  steady-state 
applications,  the  temporal  accuracy  is  not  important.  One  can  use  a  spatial  first-  , 
order  implicit  operator  (that  is,  by  setting  the  limiters  equal  zero);  i.e.,  by 
redefining  the  n's  as  !  ;  , 

«l+vum  =  {R  diag[-*(a')]  R-*}j+HAj+M      ' '  ^  V 

-    nj*+M4  =  {R  diag[-*(a*)]  R-*}k^i4Ak^j4         :     v  "^  ':  "    ^.  (3.35) 

"iw+M  =  {R  cliag[-*(a')]  R-*},+hA,^„  :  " 

Since  the  temporal  accuracy  is  not  important,  the  computation  can  be  further  ; 
reduced  by  simplifying  the  n's  diagonal  matrices,  i.e. 

"j+vw  =  {diag[-max*(a')]}j+iiAj+^4 

"j,k+vM  =  {diag[-max*(a')]}k^^4A,,+^4  ;^  (3.36) 

"ik4+M  =  {diag[-max*(a')]},+i4A,^^ 

.  If  one  applies  the  approximate  factorization  procedure  discussed  in  Section 
3.1.2  to  Eq.(3.34),  the  ADI  form  is  obtained,  i.e. 

-  *  1 

=  -M(Ej,Hw-EjW  +  (%*H,i-%-v44)  +  (Gjw+v4-Gjw.H)]"  (3.37) 


With  the  use  of  Eq.(3.33)  and  Eq.(3.34a),  the  above  equation  can  be  expressed 
in  deha  form  as  follows 

[I + Ae5  jA+ Dj]"[I + AGS    + D      +xes^C+ ]"Aq" 

=  -a(5jE  +  6„F+6^G)"  +  De  (3.38) 
where 

=  -^e(nik+wa-nj,k.M4)  (3.38a) 
=  -Ae(nj^+^-nj.,kj.v4) 

3.3.7  Application  to  thin-layer  Navier-Stokes  equations 

To  apply  the  TVD  scheme  to  the  thin-layer  Navier-Stokes  equations, 
Eq.(2.21),  is  similar  to  that  to  the  Euler  equations  [38-39].  The  viscous  term  that 
not  appears  in  the  Euler  equations  is  first  central  differenced  and  Unearlized  as 
discussed  in  Section  3.1.1  and  then  added  to  Eq.(3.38).  The  resulting  finite- 
difference  equation  is 

[I + Ae«      D^]"[I + Aes^B + D      + Ae(«  ^  C-Re'^*  ^.(J-^MJ)) + D^]"Aq" 
=  -x(6^E+s^F+s^G  -  Re-^5j.S)"+DE  (3.39) 
It  is  noted  here  that  A=At  since  A^=Aj7=Af  =  1,  and  e  is  a  parameter  in 
the  range  of  zero  to  one;  Euler  explicit  or  implicit  scheme  is  obtained  by  setting 
9=0  or  1;  if  0  =  1/2,  it  is  a  second-order  impHcit  scheme.  However,  the  order  of 
accuracy  is  degraded  due  to  the  simpUfications  that  have  been  made  in  the 
derivation  of  the  scheme.  It  is  also  noted  that  if  e  =  l/2,  Eq.(3.39)  is  similar  to 


the  Beam- Warming  scheme,  Eq.(3.10),  except  that  the  nmnerical  dissipation  tenns 
are  replaced  by  the  ones  derived  from  TVD  constraints.  In  fact,  the  TVD 
scheme  can  be  viewed  as  a  central  difference  scheme  with  a  "smart"  numerical 
dissipation  which  is  an  automatic  feed  back  mechanism  used  to  control  the 
amount  of  numerical  dissipation  [38]. 

3.3  Flow  Solver 

AppUcations  of  the  original  version  of  ARC3D  and  its  axisymmetric 
version  [40]  for  the  simulation  of  transonic  mrbulent  flow  past  a  projectile  model 
showed  that  the  codes  implemented  with  the  Baldwin-Lomax  turbulence  model 
can  give  surface  pressure  distribution  in  excellent  agreement  with  measured  data 
if  good  grid  networks  are  provided  to  the  Navier-Stokes  codes  for  the  simulation 
[25],  Nmnerical  experiments  also  showed  that  the  Beam-Warming  solution 
algorithm  is  very  sensitive  to  grid  characteristics  in  complex  flow  regions.  In 
order  to  improve  the  robustness  of  the  codes,  a  symmetric  TVD  scheme  was  first 
extended  for  modifying  the  axisymmetric  code  and  investigated  for  the  projectile 
aerodynamics  simulation;  as  reported  in  [25,41],  the  modified  code  is  indeed  less 
sensitive  to  the  grid  characteristics  in  complex  flow  regions.  Later  on,  the 
original  ARC3D  was  modified  to  implement  the  symmetric  TVD  scheme.  As 
have  been  mentioned  before,  from  the  viewpoint  of  numerical  implementation, 
the  only  difference  between  the  Beam- Warming  scheme  and  the  TVD  scheme  are 
the  numerical  dissipations.  Hence,  to  implement  the  TVD  scheme  into  the 
original  ARC3D,  one  simply  replaces  the  dissipation  operator,  Eq.(3.8),  and  the 
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dissipation  terms,  Eq,(3.9),  with  those  designed  for  the  TVD  scheme,  i.e., 
Eqs.(3.38a,  3.38b).  However,  the  latter  is  more  complicated  than  the  former. 
New  subroutines  for  computing  the  TVD  dissipations  had  been  written  and  the 
original  ARC3D  was  then  modified  into  a  TVD  thin-layer  Navier-Stokes  code 
[26].  This  TVD  code  has  been  applied  to  the  projectile  flow  problem,  and  the 
results  showed  that  the  code  was  indeed  accurate  and  robust.  Hence,  it  is 
adopted  in  this  study  for  developing  the  multi-block  computational  method. 


CHAPTER  IV 
THE  FLOW  PROBLEM  AND  SOLUTION  STRATEGY 


A  transonic  turbulent  flow  of  Mach  0.8  past  a  wing-fuselage  configuration 
shown  in  Figure  4.1,  is  adopted  as  the  model  problem  to  develop  the  multi-block 
computational  method.  The  flow  problem  considered  is  indeed  a  complex  one, 
involving  many  viscous  phenomena,  such  as  viscous  sublayer,  shock  boundary 
layer  interaction,  flow  separation,  as  well  as  complicated  boundary  geometry. 
There  are  wing  surface  pressure  measurements  available  for  assessing  the 
accuracy  of  the  numerical  solutions.  The  flow  field  is  assumed  to  be  symmetric 
with  respect  to  the  symmetry  plane  of  the  wing-fuselage  configuration  and  hence 
only  half  of  the  flow  domain  is  required  for  the  simulation,    v  , 

4.1  The  Flow  Problem  and  Preliminary  Work 

>  ■     ...  ■ . 

•  ••  -j:        '       '■  ■  ■ 
The  wing-fuselage  configuration  shown  in  Figure  4.1  is  a  model  of  a  next- 
generation  150-passenger  transport.  The  wing,  optimized  for  Mach  0.77,  has  an 
aspect  ratio  of  5.17  based  on  the  averaged  chord  length,  wing  leading  edge  sweep 
of  30",  and  an  advanced  supercritical  airfoil.  This  wing-fuselage  configuration  is 
used  as  a  model  to  develop  the  computational  methods  for  propfan  powered 
aircraft  [42].  A  rather  coarse  one  block  H-grid  for  the  configuration  was 
provided  by  the  NASA  Ames  Research  Center.  The  flow  field  is  assumed 
symmetric  with  respect  to  the  z  =  0  plane  and  hence  only  the  z  >  0  half-space  is 

45 


Figure  4.1  Wing-fuselage  configuration 
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considered  for  the  simulation.  Figure  4.2  shows  an  overview  of  the  physical 
domain  and  its  computational  domain. 

In  the  following,  the  integers  j,  k  and  1  are  used  to  denote  the  grid  points 
in  the  ^,     and  r  directions,  respectively.  Here,  ^  is  along  the  streamwise 
direction,    is  along  the  wing  spanwise  direction,  and  f  is  about  normal  to  the 
solid  surface.  The  dimensions  of  this  original  grid  network  are  jmax  =  101,  kmax 
=  45,  and  Imax  =  29.  The  solid  smfaces  are  transformed  into  part  of  the  1  =  1 
plane  with  the  shape  H,  as  shown  in  Figure  4.3;  k=l  and  k=kmax  are  symmetric 
planes;  while  j  =  l,  j=jmax,  and  l=lmax  are  artificial  outer  bovmdaries. 

To  further  illustrate  the  configuration  quantitatively,  one  takes  the 
maximum  diameter  of  the  fuselage,  D,  as  a  reference  length.  The  fuselage  is 
about  7D  long;  the  wing  span  is  about  3.3D;  the  chord  length  at  the  wing  root  is 
about  1.4D;  the  diameter  of  the  sting  is  about  0.075D;  the  distance  of  the  outer 
boimdaries  far  upstream  and  downstream  are  about  15D,  while  that  of  the  lateral 
side  is  about  lOD.  These  relative  lengths  and  the  j,  k,  1  triplet  in  physical  space 
are  shown  in  Figure  4.4.  Note  also  that  the  first  spacing  next  to  the  soUd  surface 
is  in  the  range  of  0.008D  to  0.026D. 

The  101x45x29  H-grid  one-block  system  was  first  used  for  a  preUminary 
study.  Because  of  the  specific  grid  topology,  the  thin-layer  Navier-Stokes  code 
simply  can  not  be  apphed  directly  without  modification  for  turbulent  flow 
simulation.  As  can  be  seen  from  Figure  4.3,  only  part  of  the  1  =  1  surface  is  the 
sohd  smface.  Extensive  modifications  were  made  for  proper  Jacobians  and 
metrics  calculations,  boundary  conditions,  and  eddy  viscosity  computations.  Also, 
30  grid  points  that  extremely  close  to  the  sohd  surface  were  added  to  the  original 
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Figure  42  An  overview  of  the  original  grid 
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(a)  Physical  space 
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(b)  Computational  space 
Figure  4.3  The  1=1  surface  of  the  original  grid 


Figure  4.4  Relative  lengths  in  physical  space 
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grid  to  resolve  the  thin  viscous-layer;  and  two  mirror  image  surfaces  were  also 
added  for  ease  in  symmetric  boimdary  condition  treatment.  The  101x47x59  H- 
grid  network  was  then  provided  to  the  modified  thin-layer  Navier-Stokes  code  for 
the  flow  simulation.  As  shown  in  Figure  4.5,  the  computed  results  are  in  poor 
agreement  with  experimental  data.  This  seems  to  indicate  that  the  number  of 
grid  points  on  the  wing  surface,  18  in  the  spanwise  direction  and  40  in  the 
chordwise  direction,  is  not  fine  enough  to  resolve  the  complex  flow  phenomna. 
Also,  modifications  on  the  thin-layer  Navier-Stokes  code  may  be  needed  before 
the  one-block  H-grid  topology  can  be  used  for  this  complex  turbulent  flow 
simulation.  A  flexible  and  efficient  multi-block  technique  is  more  desirable  to 
numerically  simulate  a  complex  geometry  turbulent  flow  problem,  such  as  the 
wing-fuselage  aerodynamics.  It  is  also  noted  here  that  the  calculations  of  this 
one-block  approach  were  performed  with  the  Cray  X-MP/48  of  the  NASA  Ames 
Research  Center. 

4.2  Solution  Strategy 

For  flow  past  a  complex  configuration,  such  as  the  wing-fuselage 
combination,  to  generate  an  acceptable  single  grid  for  the  entire  flow  field  is 
rather  difficult.  It  is  much  easier  to  divide  the  flow  field  into  several  subregions 
and  then  generate  a  good  grid  for  each  subregion.  Hence,  it  is  natural  to  develop 
a  solution  method  based  on  the  multi-block  grid  network  by  treating  each  block 
as  an  independent  flow  problem.  Each  block  is  then  solved  by  the  same  thin- 
layer  Navier-Stokes  code  with  the  information  between  blocks  transferred  via  the 
interface  boundary  conditions. 


Figure  4.5  Computed  results  of  the  one-block  approach 
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For  the  wing-fuselage  flow  problem,  an  effective  and  yet  simple  block-by- 
block  computational  method  is  proposed.  In  this  method,  the  entire  flow  field  is 
divided  into  six  contiguous  blocks  and  each  block  is  partially  bounded  by  a  solid 
smface.  Each  block  is  then  treated  as  an  independent  flow  problem,  which  is 
solved  by  the  same  thin-layer  Navier-Stokes  code,  with  the  interface  boundary 
conditions  updated  at  every  time  step.  This  method  does  not  need  to  worry 
about  patching  solutions  between  different  equation  zone  as  in  the  zonal  equation 
approach,  and  it  is  easier  to  deal  with  the  interface  bovmdary  conditions  than  that 
of  the  Chimera  overset  grid  approach.  However,  in  the  same  manner  as  the 
zonal  methods,  poor  interfacing  can  introduce  errors  at  the  interface  boundaries 
and  consequently  propagated  into  the  solution  field.  The  error  introduced  can 
depend  on  the  solution  behavior  and  the  grid  characteristics  around  the  interface 
as  well  as  on  the  interpolation  technique  employed  for  updating  the  interface 
boundary  conditions.  For  a  successful  multi-block  computational  method,  it  is 
then  important  and  critical  to  minimize  the  associated  interface  error  and  its 
effect  on  the  solution  accuracy. 

To  attack  the  wing-fuselage  flow  problem  by  using  the  proposed  method,  a 
flow  simulation  package  has  been  developed  which  includes  foiu-  independent 
computer  codes;  a  surface  grid  generation  (SG3D);  a  volume  grid  generation 
(GG3D);  a  multi-block  flow  solver  (WF6B);  and  a  plotting  program  (PP3D). 
They  will  be  briefly  discussed  in  the  subsequent  chapter;  however,  a  detailed 
description  can  be  found  in  a  reference  manual  [43].  Most  of  the  figures 
presented  in  this  dissertation  were  produced  by  the  plotting  program  PP3D. 
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4.3  Six-Block  Grid  Topology 

For  a  multi-block  computational  method,  the  definition  of  a  good 
composite  grid  topology  for  a  complex  configuration  turbulent  flow  simulation  is 
more  involved  and  difficult  than  that  for  an  inviscid  flow  or  a  laminar  viscous 
flow  simulation.  For  instance,  the  solid  surface  has  to  be  transformed  onto  an 
entire  plane  in  computational  space  so  that  the  thin-layer  approximation  can  be 
applied  easily;  the  Baldwin-Lomax  turbulence  model  requires  the  normal  distance 
of  each  grid  point  to  the  soUd  surface;  and  the  grids  should  be  able  to  cluster 
points  close  to  the  solid  surface  to  resolve  the  thin  viscous-layer.  Hence,  for  ease 
in  eddy  viscosity  calculation  and  in  applying  the  thin-layer  approximation,  a 
special  6-block  grid  topology  was  defined  and  investigated  for  the  wing-fuselage 
configuration.  The  main  idea  of  dividing  the  flow  field  is  to  isolate  the  wing  by  a 
cone-like  boundary  surface  that  sits  between  the  wing  and  the  fuselage,  as  shown 
in  Figure  4.6.  The  isolated  wing  was  divided  into  2  blocks  and  the  rest  of  the 
flow  field  was  divided  into  4  other  blocks.  By  doing  so,  the  flow  field  is  divided 
into  six  contiguous  blocks  such  that  each  block  is  partially  bounded  by  a  solid 
surface,  and  the  solid  surface  of  each  block  is  transformed  onto  f  =  1  plane  in 
the  computational  space.  An  overview  of  the  topology  is  shown  in  Figure  4.7 
while  Figure  4.8  gives  side  and  cross-sectional  views.  The  physical  and 
computational  space  of  each  block  is  shown  in  Figure  4.9. 

The  selected  grid  topology  has  the  clear  advantage  that  each  block  can  be 
solved  by  the  same  thin-layer  Navier-Stokes  code;  however,  the  choice  has  its 
deficiences  in  the  block  grid  generations.  For  instance,  for  the  wing  blocks.  Block 
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Figure  4.6  The  dividing  surface  isolating  the  wing 


Figure  4.7  Six-block  grid  topology  (overview) 
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Figure  4.8  Side  and  aoss-sectional  views 
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(a)  Block  1 


(b)  Block  2 

Figure  4.9  The  transformation  of  the  six  block  grids 


(d)  Block  4 
Figure  4.9  (  continued  ) 


(f)  Block  6 
Figure  4.9  (  continued  ) 
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3  and  Block  4,  shown  in  Figure  4.9(c,d),  four  of  the  six  bounding  surfaces  in  the 
computational  domain,  CLRQ,  LNTR,  RTSQ,  and  QSHC  form  a  nearly  flat 
surface  in  the  physical  space.  Therefore,  special  measure  and  care  in  grid 
generation  are  required  to  overcome  the  difficulty  resulting  from  the  drastic 
change  in  the  domain  transformation.  Also,  for  the  nose  block.  Block  1,  shown  in 
Figure  4.9(a),  the  boundary  plane  AEJF  in  the  computational  domain 
degenerated  into  a  single  line  in  physical  space  and  care  must  be  taken  when 
dealing  with  the  boundary  conditions. 


CHAPTER  V 
GRID  GENERATION 


For  the  wing-fuselage  flow  problem,  the  flow  domain  has  been  divided  into 
six  blocks,  however,  it  is  still  a  great  challenge  to  generate  acceptable  grids  for 
each  block.  To  generate  grid  networks  for  the  six-block  wing-fuselage  flow 
problem,  surface  grids  must  be  obtained  first.  A  special  surface  grid  generation 
code  (SG3D)  was  developed  to  utilize  the  original  H-grid  network  for 
constructing  the  required  surface  grids.  The  resulting  surface  grids  will  be 
presented  in  section  5.1.    Generating  volume  grids  for  the  six-block  flow  domain, 
it  is  required  to  develop  an  adequate  grid  generation  code.  Various  grid 
generation  methods  [44-51]  have  been  studied  in  that  regards.  The  main  concern 
is  that  the  code  can  produce  grid  networks  that  satisfy  the  required  grid 
characteristics,  such  as:  (a)  grid  lines  in  different  directions  must  be  nearly 
orthogonal  to  each  other;  (b)  grid  lines  in  {--direction  must  be  about  normal  to 
the  solid  surface;  (c)  r  =  constant  grid  surfaces  must  be  clustered  near  the  solid 
surface;  (d)  grid  lines  must  be  continuous  across  the  interfaces;  and  (e)  the  blocks 
must  be  contiguous  (no  holes  or  overlapping).  A  combination  of  two-boundary 
algebraic  and  elliptic  grid  generation  methods  will  satisfy  these  requirements  best. 
Accordingly,  a  general  purpose  3-D  Grid  Generation  code  (GG3D)  that  based  on 
these  two  methods  was  developed  and  applied  to  generate  the  3-D  grid  networks. 
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A  detailed  description  of  GG3D  can  be  found  in  [43].  The  techniques  employed 
to  generate  the  six  block  grids  are  described  in  Section  5.2  and  the  resulting  grid 
networks  are  presented  in  Section  5.3. 

5.1  Surface  Grid  Generation 

To  generate  grid  networks  for  the  six  block  grids,  surface  grids  must  be 
obtained  first.  A  two-boundary  algebraic  grid  generation  method  is  used  to 
generate  volume  grids;  consequently,  only  the  solid  surface  grid  (1=1)  and  the 
cooresponding  outer  surface  grid  (l=lmax)  are  required.  The  other  boundary 
surfaces  are  then  defined  by  the  straight  Unes  that  connect  the  boundary  curves  of 
the  solid  and  outer  surface  grids. 

A  special  surface  grid  generation  code  (SG3D)  has  been  developed  that 
uses  the  original  H-grid  to  construct  each  surface  grids  by  cubic  splining, 
transfinite  interpolation,  strentching,  shifting  and  dividing.    The  resulting  surface 
grids  obtained  are  shown  in  Figure  5.1;  the  dimension  of  these  surface  grids  are 
28x23,  40x12,  40x35,  40x35,  40x12  and  26x23  for  Block  1-6,  respectively. 

It  is  noted  here  that  the  outer  surfaces  are  defined  according  to  the  grid 
topology  and  will  have  strong  effect  on  the  grid  characteristics  of  the  block  grid  to 
be  generated.  Hence,  the  boundary  curves  of  each  outer  surface  should  be 
defined  carefully.  For  block  1,  a  cone-like  shape  that  is  resemble  to  the  solid 
surface  is  desired  for  the  outer  surface.  As  shown  in  the  grid  topology.  Figure 
4.7,  the  boundary  curve  GHI  is  defined  such  that  the  bounding  surface  lies  in  the 
middle  of  the  wing  leading  edge  and  the  fore-body  of  the  fuselage.  Similarly,  the 
outer  surface  of  block  6  is  defined  such  that  bounding  surface  of  the  boundary 


(a)  Solid  surfaces 
Figure  5.1  The  6-bIock  surface  grids 
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curve  MNO  lies  in  the  middle  of  the  wing  trailing  edge  and  the  aft-body.  The 
boundary  curves  Hb'N  of  block  2  and  Ha'N  of  block  5  are  also  defined  such  that 
the  boundary  surfaces  Hb'NLC  and  Ha'NLC  lie  in  the  middle  of  the  wing  and  the 
fuselage.  For  the  outer  surfaces  of  block  3  and  block  4,  the  boundary  curve  HN, 
that  is  defined  as  a  straight  line,  is  located  such  that  the  interface  surface  of  block 
3  and  block  4,  i.e.,  surface  HRCQLN,  is  not  incline  to  either  upper  side  or  lower 
side  of  the  wing.  In  other  words,  the  plane  spanned  by  HCLN  should  pass 
throught  the  wing  tip,  and  hence  divides  the  wing  into  upper  wing  and  lower  ; 
wing.  If  this  were  not  the  case,  one  would  have  difficulty  to  generate  an  usable 
volume  grid  for  either  block  3  or  block  4. 

»   '  '■  . 

5.2  Volume  Grid  Generation  ■ 

■  'C 

*  '.       -.  -"  ^ 

5.2.1  Two-Boundary  Grid  Generation 

The  two-boundary  grid  generation  technique  used  is  based  on  Hermite 
transfinite  interpolation  [46-48].  Two  surface  grids  are  specified  and  the  interior 
grid  lines  are  connected  by  curves  such  that  boundary  orthogonality  is  enforced. 
The  other  four  boundary  surfaces  are  then  defined  using  straight  lines.  The 
distribution  of  interior  stufaces  is  defined  by  a  suitable  clustering  function  such 
that  more  grid  points  are  clustered  near  the  solid  surface  as  it  is  usually  required. 
A  description  of  clustering  functions  is  given  in  Appendix  C. 

Denoting  the  position  vectors  of  the  two  specified  boundary  siufaces  as 
r(?,r;,0)  =  r,(^,,,)    and      r(e,»?,l)  =  r2(^,,,)  (5.1) 
and  the  derivatives  with  respective  to  f  as 
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ar 

the  following  functional  form,  based  on  Hermite  transfinite  interpolation,  may  be 
written  for  the  connecting  curves: 

ra.'/.r)  =  ri(|,„)hi(f)  +  r2(^,„)h2(f)  +  DRih3(r)  +DR2h4(f)  (5.3) 
To  satisfy  Eq.(5.1)  and  Eq.(5.2),  the  following  must  be  true: 

h.(0)  =  1,  h.(i)  =  0,  iMOL  =  0,  iMlI  =  0 

h,(0)  =  0.  h,(l)  .1.^=0,    ^  =  0 

h,(o)  =  0,  h3(i)  =  0,  ^M2L  =  1,  iMlL  =  0 

^  ar  ar 

h,(0)  =  0,  h,(i)  ,0,^=0,  iMl).  =  1 

af  ar 

Hence,  a  cubic  polynomial  in  {•  for  hj(f )  satisfying  conditions  given  in  Eq,(5.4) 
can  be  determined.  They  are 

hi  =  2r^-3r^+l 

h2  =  -2r'+3f2 

ha  =  r'-2r'+r  (5.5) 

Here,  f  is  a  normalized  point  distribution,  and  is  usually  defined  by  a  clustering 
function  as  shown  in  Appendix  C. 

The  expressions  for  the  derivatives,  DR,  and  DRj,  appearing  in  Eq.(5.3) 
are  the  key  elements  used  to  enforce  the  boundary  orthogonahty  and  are  derived 
as  follows. 
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A  vector  normal  to  the  surface  r  =  0  that  points  into  the  spatial  domain  is 
given  by 

where  x  denotes  the  vector  cross  product.  A  vector  that  tangents  to  the 
connecting  curve  at  r=0  can  be  defined  as 

ar(e,>?.0) 
e  =   

ac 

Hence,  in  order  for  the  connecting  curves  to  be  normal  to  the  r  =  0  surface,  the 
cross  product  of  the  vector  tangent  to  the  connecting  curves  e  and  the  surface 
normal  n  must  be  zero,  i.e., 

exn  =  Oatf=0  (5.5) 
or  written  in  component  form 

r  .^ci^_-^  \y 
ar    a»7   ai    ac   ari      ar    at;  a^  "aT  ^ 

af  ^  ar;  ac  '  af  at;  ^"  ar    a»;  a^  "ai~  "a^^ 

+  f  iL.  ( IL.  .££_        ^  :?y_  /  iy_ 

ar    ar?  a^    ac  ar?     ar    ar?  ae  . 
=  0  (5.5a) 
From  the  above  equation,  it  can  be  seen  that  the  connecting  curves  will  be 
perpendicular  to  the  surface  r  =  0  if  the  following  conditions  are  satisfied  at  f  = 
0: 

DX,  =  ^  =  -CK,(e,.)(  JLJL.^JL^ 
ar  ^^^"^^  dr,  ar,  ' 
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DYi  -  —  -  CKi(e,»7)(-— — --— — )  -    '  (5.7) 

ari  di  arf 

-        -  rT<r      \(        3y    ax  ay 

Here,  CK^  are  arbitrary  functions  introduced  to  control  the  shape  of  the 
connecting  curves,  and  if  set  to  zero  the  curves  become  straight  lines. 

In  a  similar  manner,  one  can  show  that  the  connecting  curves  will  intersect 
the  f  =  1  surface  perpendicularly  if 

af  drj  dE,  or) 

DY,  -  —  -  CK,(c,.)(  ——-—_)  (5.8) 

ny  _        _  ax   ay    ax  ay 

DZ,  -  —  - -CK,(|,,)(— — -— — ) 

are  satisfied  at  r  =  1. 

By  substituting  Eqs.(5.7-5.8)  into  Eq.(5.3),  one  obtains  the  following 
expressions  for  the  resulting  volume  grid  with  boundary  orthogonality  enforced: 

X(^,»7,f)  =  Xi(^,„)hi(r)  +  X2(^,r7)h2(r)  +  DXih3(0  +  DX2h4(r) 

y(?,'/,f)  =  yi(^,'7)hi(f)+y2a,r7)h2(f)+DYih3(f)+DY2h4(r)       -  (5.9) 

zU.'y.f)  =  Zi(?,'7)hi(f)  +  Z2(e,r/)h2(f)  +  DZ,h3(r)  +  DZ2h4(f)  ' 

It  is  noted  here  that  the  derivatives  on  the  right  hand  side  of  Eqs.(5.7-5.8)  are 

discretized  by  proper  finite-  difference  fomula,  e.g.,  central  difference 

3yi  _  yi(k+l)-yi(k-l)  ,  '  ' 

a,   -  .(k-fl).,(k-l)     ^  ^(^^  > 

and  C  and  ri  are  normalized  as  follows: 
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CG)  =  (j-l)/(jmax-l),    j  =  l,2,...,jmax 

r,(k)  =  (k-l)/(kmax-l),  k  =  l,2,...,kmax  (5.10) 
The  CK's  appearing  in  Eqs.(5.7-5.8)  have  a  strong  effect  on  the  shape  of 
the  connecting  curves.  In  fact,  they  control  the  depth  of  the  perpendicular  part  of 
the  curves.  Values  for  the  CK's  are  usually  specified  as  constants  by  trial  and 
error.  A  shape  function  was  devised  to  modify  the  CK's  so  that  different  values 
can  be  specified  at  different  j,  k  locations.  The  shape  function  is  of  the  form 
CK  =  CKI*cosx{CPJ[0-l)/Omax-l)]  +  CSJ}  > 

«cos^{CPK[(k-l)/(kmax-l)]  +  CSK}  (5.11) 
where  CKI  is  a  constant  to  be  specified,  CPJ  and  CPK  are  used  to  control  the 
frequency  while  CSJ  and  CSK  are  used  to  control  the  phase  of  the  relative  cosine 
function.  In  doing  so,  values  for  the  CK's  become  functions  of  j  and  k,  and  can 
be  smoothly  distributed  on  the  surface.  For  example,  if  one  specifies  CPJ  =  0.5, 
CSJ  =  0,  then  CK=CKI  when  j  =  l  and  CK  decreases  as  j  increases;  or  if  one 
specifies  CPJ  =  1,  CSJ =-0.5  then  CK  increases  from  zero  to  CKI  and  then 
decreases  to  zero  as  j  varies  fi-om  1  to  jmax.  Figure  5.2  shows  some  examples  to 
illustrate  this  shape  function. 

5.2.2  Elliptic  Grid  Generation 

The  technique  and  subroutines  employed  are  based  on  the  elliptic  grid 
generation  system  of  Thompson  [49-51].  The  volume  grid  is  constructed  by 
forcing  the  Cartesian  coordinates  of  interior  points  to  satisfy  specially  devised 
Poisson  equations.  In  the  Poisson  equations,  three  control  functions  for  each 
curvilinear  direction  are  set  up  to  control  the  grid  point  distributions.  The  system 
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CP  j-e.  5 ,  cs  j-e ,  cPK-e ,  csK-e        cp  >e ,  cs  j-0 .  cpk-  .  s ,  csk— .  4 


CPU-.  8 ,  cs  J— a .  4,  CPK- .  5,  csK— .  5  CP  j-e,  cs  j-e,  cPK-a,  csK-9 


Figure  5.2  Shape  function  for  CK's 
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of  Poisson  equations  are  solved  by  point  SOR  (success  overrelaxation)  iterations. 
The  user-provided  initial  volume  grid  is  used  not  only  to  start  the  SOR  iterations 
but  also  to  set  up  the  control  functions. 

The  elliptic  grid  generation  system  is  based  on  the  Poisson  equations  of 
the  form 

=  (v^.v^Pi 

v^r,  =  (v„.v„)P2  (5.13) 

A  =  (vr-vr)P3 

where  the  P's  are  so-called  control  functions  which  can  be  specified  to  control  the 
grid  smoothness  and/or  orthogonahty.  Since  the  actual  computations  are 
performed  in  the  computational  domain  where  ^,  rt  and  $•  are  the  independent 
variables,  with  r  =  (x,y,z)^  as  the  dependent  variables.  The  poisson  equations, 
Eq.(5.13),  are  then  transformed  to  computational  domain.  The  transformed 
equations  £u-e 


a^r  a^r  d^v 

(V^.V^)  __  +  (V„.V„)  —  +  (Vr.7„)  -—  -H 
a^OT)  dtjdr}  d^dr/ 

a^r  a^r  a^r 


acsr  a»7a$-  '  a^ar 

(7e.vOPi^+(Vr,.V„)P2^+(Vf.Vf)P3^  =  0  (5.14) 
OS  or/  af 

The  finite-difference  equafions  of  Eq.(5.14)  are  obtained  by  replacing  the 
derivatives  with  central  differences  and  these  equations  are  then  solved  by  point 
SOR  iteration  procedure. 
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If  the  control  function  P's  are  evaluated  from  the  initial  grid,  the  three 
components  of  the  elliptic  grid  generation  system,  Eq.(5.14),  provide  a  set  of 
three  equations  in  which  the  control  functions  P's  are  treated  as  unknowns  and 
the  values  for  the  position  vectors  are  obtained  from  the  initial  grid.  The 
transformed  poisson  equations,  Eq.(5.14),  can  be  rewritten  as 

(ve- vOPi +  (v»7  •  v,,)P2  ^ + (vr .  vr)P3 


d^d^  dfidi  arac 

/       X  /       X  ,       ,  aV 

(VC-Vr,)—  +(V„.V„)— —  +(vr.V„)—  + 
a^drj  drjdri  d^drj 

which  can  be  solved  simultaneously  at  each  point  for  the  three  control  functions, 
Pj,  i=l,3.  If  the  control  functions  are  computed  solely  based  on  the  above 
equations,  and  the  elliptic  system,  Eq.(5.14),  will  reproduce  the  same  grid  as  the 
initial  grid  which  is  of  trivial  interest.  However,  if  the  computed  control  functions 
are  smoothed  before  they  are  applied  into  the  elliptic  system,  a  smoother  grid 
network  can  then  be  produced. 

A  variation  to  the  evaluation  of  the  control  functions  is  obtained  by 
assuming  that  the  grid  lines  are  orthogonal  to  each  other,  and  Eq.(5.15)  is 
reduced  to 

(ve.vOPi-^+(V,.V,)P2i^^+(Vr.Vf)P3^  = 

oT)  a$" 

-[("^^'^O  4hc  +  (^''•^'')  -fr  +  (^^-"^f )  -TT-  ]  (5.16) 

aca?  dr}dri  acaj-  ^ 


The  control  functions  evaluated  from  these  equations  will  force  the  elliptic  system 
to  produce  a  more  orthogonal  grid. 

5.3  The  Six-Block  Grid 

A  set  of  volume  grids  with  60  grid  points  in  the  I-direction  was  first 
generated  by  two-boundary  grid  generation  with  grid  hues  set  normal  to  the  solid 
surface  whenever  possible.  The  first  spacing  off  of  the  solid  surface  was  set  to 
about  2xlO'^D,  which  is  considered  fine  enough  for  turbulence  simulations. 
Figure  5.3  shows  the  3-D  close-up  plots  of  the  six  block  grids  with  cut-offs  to 
show  the  details  of  the  interior  surfaces.  These  grids  were  generated  from  the 
soUd  and  outer  surfaces  described  in  Section  5.1.  The  other  four  boundary 
surfaces  were  then  defined  by  straight  lines  that  connect  the  boundary  curves  of 
the  solid  and  outer  surfaces.  Sixty  grid  points  were  distributed  on  each  of  the 
connecting  straight  lines  by  using  the  hyperboHc  tangent  clustering  function  with 
the  first  spacing  of  2xlO'*D  and  the  last  spacing  of  six  times  of  the  average 
spacing,  and  hence,  more  grid  points  were  clustered  near  the  solid  surfaces.  To 
generate  grid  network  using  two-boundary  grid  generation,  one  has  to  face  a 
difficulty  in  specifying  a  proper  value  for  the  parameter  CK,  Eqs.(5.7-5.8).  In 
general,  the  value  for  CK  is  specified  by  trial  and  error,  however,  its  variation  on 
a  surface  can  be  justified  by  the  use  of  the  shape  function,  Eq.(5.11),  according  to 
the  geometry  of  the  grid  to  be  generated.  For  the  six  block  grids,  the  grid 
orthogonahty  on  the  outer  surface  was  not  enforced,  that  is  CK2  =  0,  to  avoid  the 
possibiUty  of  grid  overiapping.  On  the  other  hand,  the  grid  orthogonality  on  the 


(b)  Block  2  -  Figure  shows  the  surfaces 

1=1, 1=50,  k=3,  k=13,  j  =  17 


Figure  53  Six  block  grids  - 


-  3D  close-up  with  cut-off 
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(c)  Block  3  ~  Figure  shows  the  surfaces 
1=1,  1=50,  k=l,  k=15,  k=30,  j  =  10,  j  =  18,  j=32 


Figure  5.3  (continued) 
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(e)  Block  5  --  Figure  shows  the  surfaces 

1=1,  1=50,  k=3,  k=13,  j  =  17 


(f)  Block  6  ~  Figure  shows  the  surfaces 

1=1,  1=50,  k=5,k=15,j=l,j=10 


Figure  53  (continued) 
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solid  surface  is  desirable.  From  numerical  experiments  conducted  the  following 
was  one  of  the  best  to  set  up  the  CKj  distribution  by  the  shape  function  for  each 
block. 


Block  1  - 

CPJ 

=0.5 

CSJ= 

0 

CPK= 

0 

CSK= 

0 

Block  2  - 

CPJ 

=  1 

CSJ  = 

-0.5 

CPK= 

0.5 

CSK= 

0 

Block  3  - 

CPJ 

=  0.8 

CSJ  = 

-0.4 

CPK= 

1 

CSK= 

-0.5 

Block  4  - 

CPJ 

=0.8 

CSJ  = 

-0.4 

CPK= 

1 

CSK= 

-0.5 

Block  5  - 

CPJ 

=  1 

CSJ  = 

-0.5 

CPK= 

0.5 

CSK= 

0 

Block  6  - 

CPJ 

=  0.5 

CSJ  = 

-0.5 

CPK= 

0 

CSK= 

0 

These  distributions  are  shown  in  Figure  5.4.  For  block  1,  CKj  decreases  as  j 
increases  and  becomes  zero  when  j  =  jmax.  This  is  reasonable  since  the  j  =  jmax 
surface  incline  about  45  degrees  to  the  solid  surface.  Similar  argument  can  be 
addressed  for  the  other  blocks.  To  further  improve  the  smoothness  of  the  grid 
network,  these  algebraic  grids  were  then  used  as  initial  grids  for  elliptic  grid 
generation.  Smoother  grids  were  obtained  for  block  1,  2,  5,  and  6.  However,  it 
failed  to  generate  elliptic  grid  for  either  block  3  or  4.  This  may  due  to  their 
special  geometries  such  as  three  boimdary  surfaces  form  nearly  a  flat  plane. 
Figure  5.5  illustrates  the  difference  between  algebraic  and  elliptic  grids.  Notice 
that  the  elliptic  grids  are  smoother  and  more  orthogonal  than  the  algebraic  grids. 

These  six  block  grids  were  further  modified  to  include  mirror  image  planes 
for  imposing  symmetric  boundary  conditions  in  the  Navier-Stokes  simulations. 
The  resulting  grid  network  was  then  used  as  a  basic  grid  in  this  stody  and  is 
referred  to  as  the  "base  grid"  hereafter.  The  dimensions  of  this  base  grid  are 
28x25x60,  40x13x60,  40x35x60,  40x35x60,  40x13x60,  and  26x25x60  for  blocks  1-6, 


Figure  5,4  Shape  function  for  the  six-block  grid 
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Block  1,  k-13,  fllgsbralc  Grid 


Z 


Block  1,  k-13.  Elliptic  Grid 


(a)  Block  1  -  k=13 
Figure  5.5  Elliptic  and  algebraic  grids 
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(b)Block6-k=13 
Figure  5.5  (continued) 
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respectively.  Figure  5.6  shows  some  detailed  plots  to  illustrate  the  grid 
characteristics.  Figures  5.6a-d  show  that  the  grid  lines  are  about  normal  to  solid 
surface  whenever  possible.  Figures  5.6e-f  show  that  some  k= constant  surfaces 
are  quite  skew  for  block  2,  3,  4,  and  5.  Figure  5.6g  shows  the  grid  distribution  at 
k  =  l,  12,  and  18  wing  cross-sections. 

In  order  to  investigate  the  effect  of  grid  resolution  on  the  flow  solution,  a 
series  of  grid  refinement  studies  have  been  performed.  Several  locally  refined 
grids  (near  the  wing  tip  or  near  the  soUd  surface,  etc.)  were  generated  and  used 
to  smdy  the  possible  effect  due  to  the  grid  characteristics.     A  set  of  "refined 
grids"  was  also  generated  by  doubling  the  grid  points  in  the  wing  core-wise 
direction.  The  dimensions  of  these  refined  grids  are  28x25x60,  79x13x60, 
79x35x60,  79x35x60,  79x13x60  and  26x25x60.  Figure  5.7  shows  the  solid  surfaces 
of  block  2  and  3  of  this  refined  grid.  Note  that  the  grid  lines  were  generated  by 
using  a  cubic  sphne  so  that  the  soUd  surfaces  would  not  be  distorted. 
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L 


j-1  or  48 


Block  3  or  4 
k-35 


J-40  or  1 


Slock  1  k-13 


Block  6  k-13 


(a)  Block  1,  3,  6 


(b)  Qose-up  to  show  the  boundary  orthogonality 
Figure  5.6  Some  detail  plots  for  six-block  grids 
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Block  2  Block  5 

J-20  j-2e 


(c)  Block  2,  3,  4,  5 


(d)  Close-up  to  show  the  boundary  orthogonality 
Figure  5.6  (continued) 
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(e)  Block  3,  4  -  k  surfaces  are  quite  skew 


(f)  Block  2,  5  -  k  surfaces  are  quite  skew 
Figure  5.6  (continued) 


86 


Slock  4  k-18 


Block  3  k-18 


Block  4  k-12 


Block  3    k-1  1-1-35 


(g)  Block  3,  4  -  k=l,  12,  18  to  show  the  grid  distribution 
on  the  leading  and  trailing  edge 


Figure  5.6  (continued) 
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Base  Grid  -  48x35  on  uing  surfaca 


(a)  Base  grid  -  Block  3,  4  -  40x35 
Block  2,  5  --  40x13 


Refined  Grid  -  79x35  on  wing  surface 


(b)  Refined  grid  -  Block  3,  4  -  79x13 
Block  2,  5  -  79x13 

Figure  5.7  Base  grid  and  refined  grid 


CHAPTER  VI 
SOLUTION  METHOD 


The  multi-block  computational  method  considered  is  rather  simple  and 
straightforward  in  application.  For  the  wing-fuselage  flow  problem,  the  flow  field 
has  been  properly  divided  into  six  contiguous  blocks,  and  the  grid  network  has 
also  been  generated.  In  the  solution  process,  each  block  is  then  treated  as  an 
independent  flow  problem  that  is  solved  by  the  same  thin-layer  Navier-Stokes 
code.  The  interface  boimdary  conditions  between  blocks  are  updated  with  a 
distance  weighted  interpolation  before  the  computations  march  for  another  time 
step.  A  computer  code  has  been  developed  and  investigated  for  the  six-block 
wing-fuselage  aerodynamics  in  which  the  unsteady  3-D  TVD  thin-layer  Navier- 
Stokes  code  is  made  to  a  subroutine  subprogram.  Since  the  multi-block  computer 
code  has  to  deal  with  several  blocks  of  different  configurations  and  their 
interfaces,  the  code  has  to  be  modular  and  flexible.  Also,  the  code  has  to  adapt 
to  the  architecture  of  the  host  computer  system.  A  six-block  wing-fuselage 
Navier-Stokes  code  (WF6B)  was  then  written  for  the  available  computer  system, 
Cray-2  of  NAS  system  at  NASA  Ames  Research  Center.  The  special  feature  of 
the  supercomputer,  Cray-2,  is  its  huge  main  memory,  however,  the  computing 
speed  is  a  Uttle  slower  than  Cray  Y-MP.  In  order  to  take  advantage  of  the  large 
main  memory  while  not  sacrijSce  the  computing  speed,  WF6B  uses  part  of  the 


88 


89 

main  memory  for  storing  the  intermediate  data.  The  other  main  memory  is  then 
used  as  active  area  for  executing  computations.  A  set  of  common  arrays  is  then 
allocated  at  the  active  area  and  used  for  the  computations  of  each  blocks.  A 
block  diagram  to  show  this  data  arrangement  is  given  in  Figure  6.1.  A  detailed 
description  of  WF6B  can  be  found  in  [43]. 

In  order  to  illustrate  the  multi-block  solution  method  a  flowchart  of  WF6B 
is  given  in  Figure  6,2.  As  shown  in  the  figure,  two  do-loops  run  over  each  blocks 
and  the  nmnber  of  time  steps  requested.  In  the  inner  loop,  for  each  block,  BC  is 
used  to  set  up  the  boundary  conditions  and  then  TLNS  is  to  formulate  the 
governing  equations  and  perform  the  three  inversions.  A  value  of  L2  residual  for 
each  block  is  computed  by  RESID  every  10  time  steps.  After  the  computations 
run  over  each  block,  INTERBC  updates  the  interface  boundary  conditions,  and 
then  the  computation  continues  to  the  next  time  step.  Finally,  the  pressure  and 
the  shear  stress  coefficients  are  computed  by  CPP  and  CSP,  respectively.  In  what 
follows,  the  techniques  used  to  calculate  the  metric  coefficients  and  Jacobians  of 
transformation  are  described  first.  Then  the  boundary  conditions  imposed  on  the 
boundary  surfaces  of  the  six-block  flow  problem  are  discussed  while  the  interface 
boundary  conditions  are  given  in  a  separate  section.  Fmally,  the  equations  used 
to  compute  the  12  residual  and  pressure  and  shear  stress  coefficients  are  given. 

6.1  Jacobian  and  Metric  Coefficient 

The  metric  coefficients  and  Jacobians  of  transformation  can  be  calculated 
from  Eqs.(2.16-2.17),  if  the  covariant  metrics,  x^,  y^,  z^,  etc,  are  known.  For  a 
second-order  approximation,  the  covariant  metrics  are  calculated  at  interior  grid 
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Block  1 

28x25x60 
rl(3, j,k,l) 
ql(j,k,l,6) 


Block  2 
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r2(3, j,k,l) 
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Block  4 
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r5(3, j,k,l) 
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results  sweep  out 
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Figure  6.1  Data  arrangement  in  WF6B 
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Figure  6.2  The  flowchart  of  WF6B 
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points  by  a  central-difference  formula  and  second-order  one-sided  differencing  at 
the  boundaries,  e.g., 

=  (Xj+ijc4-XiW/2  for  1  <  j  <  jmax 

Xf  =  (-3Xj,m+4Xj+i^-Xj+2j^/2        for  j  =  1  (6.1) 

=  (  3Xj^-4Xj.ij^+Xj.2j^/2  for  j  =  jmax 

The  Jacobian,  J,  of  each  grid  cell  is  then  computed  from  Eq.(2.17).  The 
contravariant  metric  coefficients  are  defined  in  Eq.(2.16),  and  can  be  calculated 
accordingly. 

The  Jacobians  and  the  metric  coefficients  evaluated  by  the  above  finite- 
difference  approach  are  rather  sensitive  to  grid  characteristics,  and  can  produce  a 
large  truncation  error  if  the  grid  network  does  not  have  good  properties  of 
smoothness  and  orthogonality.  It  is  known  that  they  can  also  be  approximated 
from  geometric  relations  such  as  cell  volume  and  cell-face  area.  In  the         ,  ^ 
discretized  sense,  the  Jacobian  defined  in  Eq.(2.17)  is  equivalent  to  one  over  the 
cell  volume,  1/Vol,  and  the  contravariant  metric  coefficients  defined  in  Eq.(2.16) 
are  equivalent  to  one  over  the  cell-face  normal.  If  one  defines  a  computational 
cell  that  is  centered  at  the  corresponding  grid  point,  with  its  cell-faces  located 
half-way  to  the  neighboring  points,  it  is  then  appropriate  to  replace  the  Jacobian, 
J,  by  the  corresponding  cell  density,  1/Vol.  Furthermore,  since  C^/J,  ^y/J,  and 
iji  defined  the  x,  y,  z  components  of  the  normals  (not  unit  normals)  to  the  local 
tangent  plane  of  the  constant  |  surface,  one  can  replace  by  n^/Wo\, 

Uj-y/Vol,  and  Uj^/Vol,  where  n^^  Ujy,  and      are  the  components  of  the  normals  to 
the  local  constant  j  planes.  Likewise,  one  can  replace     rj^  t)^  and  c,,  fy,  by 
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n^/Vol,  n,jy/Vol,  n,„/Vol  and  n^j/Vol,  n,y/Vol,  n,^/Vol,  respectively. 

A  general  hexahedron,  defined  by  eight  arbitrary  comer  points,  can  be 
partitioned  into  six  tetrahedrons,  Figure  6.3a.  The  volume  of  a  tetrahedron 
denoted  by  its  vertices  (a,b,c,d)  equals  one-sixth  of  the  triple  product  of  the  three 
vectors  emanating  from  one  of  the  vertices  and  ordered  according  to  the  right- 
hand  directions.  Figure  6.3b.  For  example, 

V'*(a,b,c,d)  =  [r.,.(r^  x  r,,)]/6  (6.2) 
=  [(Xc-Xd)(yb-yd)(z,-Zd)  +  (Xb-Xd)(y.-yd)(z,-Zd) 

+  (Xa-Xd)(yc-yd)(Zb-Zd)-(vxd)(yb-yd)(Zc-Zd) 
-  (vxd)(ya-yd)(Zb-Zd)-(vxd)(yc-yd)(ZaZd)]/6 

where     =  Tj  -  Fj,  and  Tj  and  Tj  are  the  position  vectors.  If  the  eight  vertices  of 
the  hexahedron  are  denoted  by  numerals  1-8  and  are  associated  with  the  j,k,l 
triplets  as  follows: 

l-j,k,l  2-j  +  l,k,l 

3-j+l,k+l,l  4-j,k+l,l 

5-j,k,l+l  6-j  +  l,k,l+l 

7  -  j+l,k+l,l+l  8-j,k+l,l+l 
the  volume  of  the  hexahedron  is 

Voiak,l)  =  V*'*(l,2,3,7)  +  V***(l,2,7,6)  +  V***(l,5,6,7) 

+  V*'*(l,3,4,7)  +  V*(l,4,8,7)  +  V*''(l,5,7,8)  (6.3) 
The  Jacobian  of  a  grid-centered  computational  cell,  which  is  bounded  by  the 
surfaces  at  the  middle  of  two  adjacent  grid  points,  is  computed  as  one  over  the 
average  volmne  of  all  the  associated  hexahedrons. 
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Figure  6.3    Volume  computation  of  a  hexahedron 
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The  metric  coefficients,  based  on  the  cell-face  area,  are  calculated  in  a 
similar  maimer.  The  area  vector  of  a  triangle  denoted  by  its  vertices  (a,b,c)  is 
A^(a,b,c)  =  (r^b  X  r^)/2 

=  {[(yb-yJ(vzb)-(yb-ya)(Zc-Zb)]» 

+  [(Zb-Za)(VXb)-(Zb-Za)(VXb)]j 

+  [(vxa)(yc-yb)-(vxa)(yc-yb)]k}/2  (6.4) 

The  area  vector  of  a  general  surface,  defined  by  four  comer  points  (a,b,c,d),  may 
be  defined  as 

AREA(a,b,c,d)  =  A^(a,b,c)  +  A''(c,d,a)  (6.5) 
Hence,  nj^  are  the  x,y,z  components  of  AREA(a,b,c,d).  For  a  constant  j 

plane,  one  has 

%x  =  {[(yb-ya)(Zc-Zb)-(yb-y.)(Zc-Zb)] 
+  [(yd-yc)(vzd)-(yd-ye)(vZd)]}/2 

njy  =  {[(Zb-Za)(Xc-Xb)-(Zb-Za)(x,-Xb)] 

+  [(Zd-Zc)(vXd)-(Zd-Zc)(vXd)]}/2  /•  (6.6) 

Djz  =  {[(vxa)(yc-ybHvxJ(yc-yb)]  -  - 

+  [(Xd-Xc)(ya-yd)-(Xd-Xc)(ya-yd)]}/2 

and  similarly,  one  can  compute  n^^  n^y,  n^,  and  n^^  n,y,  n,^. 

6.2  Boundary  Condition 

There  are  different  types  of  boundaries  encountered  in  nmnerical 
simulations,  such  as  solid  surface,  far-field,  inflow/outflow,  and  symmetric 
boundaries.  The  type  of  boundary  condition  specified  at  a  boundary  surface 
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depends  on  the  physical  behavior  and  numerical  requirement  at  the  surface. 
According  to  the  grid  topology,  the  proper  boundary  condition  has  been  set  up 
for  each  block  in  WF6B.  Figure  6.4  shows  the  boundary  conditions  applied  at  the 
boundary  surfaces  of  each  block  and  the  relative  interfaces  between  them.  In  the 
figure,  SOLIDBC,  FSBC,  SYMBC,  SGLRBC,  OUTBC  and  INTERBC  stand  for 
soUd,  freestream,  symmetric,  singular,  outflow  and  interface  boundary  condition, 
respectively.  The  j  =  1  surface  of  block  1  is  a  special  surface  that  degenerates  into 
a  single  line  in  physical  space,  special  treatment  is  required.  Since  on  the  surface, 
each  ri  line  maps  to  a  point  in  physical  space,  the  same  boundary  condition  is  set 
at  each  point  on  this  line,  and  the  value  is  obtained  by  distance  weighted 
extrapolation  from  the  interior  points.  The  other  boundary  conditions  are 
discussed  in  the  following  and  the  interface  boundary  condition  is  given  in  the 
next  section. 

6.2.1  Solid  Boundary  Condition 

In  the  computer  code,  the  r  =  0  plane  is  assumed  the  solid  surface.  The 
techniques  employed  to  obtain  the  velocities,  pressure,  density,  and  total  energy 
will  be  described  for  both  inviscid  and  viscous  flows,  and  then  a  slow  start 
mechanism  for  viscous  flows  is  given. 

Inviscid  flow.  At  a  soUd  surface,  the  flow  tangency  condition  is  satisfied. 
So  the  contravariant  velocity  component  in  f -direction  W  is  set  to  be  zero  for 
non-injecting  wall  and  the  other  contravariant  velocities  U  and  V  are  foimd  by 
extrapolation,  e.g.. 
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Vjw  =  2V^^  -  V.,^  (6.7) 

The  velocity  components  u,  v  and  w  are  then  obtained  from  Eq.(2.20c).  Note 
that  ^, ,  r;,  and    are  zero  for  a  fixed  grid  system  employed  here. 

The  total  energy  E  at  the  solid  surface  is  obtained  from  Eq.(2.7)  in  which 
the  pressure  is  computed  from  the  normal  momentum  equation.  The  normal 
momentum  equation  is  obtained  by  combining  the  three  transformed  inviscid 
momentimi  equations  [10],  i.e.. 


a(pu/J)  1 

a[(puU+^j))/J]  1 

a[(puV+r,J))/J]  _ 

a[(puW+rj5)/J] 

=  0 

at 

a(pv/J)  ^ 

3[(pvU  +  e,p)/J]  ^ 

a[(/.vV+„,p)/J]  ^ 

a[(pvW+r,p)/J] 

=  0 

at 

ae 

ar; 

ar 

a(pw/J)  1 

a[(pwU+ej))/J]  1 

a[(pwV+r,j))/J]  1 

a[(pwW+rj))/J] 

=  0 

at 

a€ 

a»? 

ar 

If  combining  r^'C-momentum  +  ry •'/-momentum  +  r^'r-momentum  with  the 
use  of  the  continuity  equation  and  metric  identities,  Eq.(2.19),  one  has  the  normal 
momentum  equation 

Uxrx+ Cyry + i.r,)Pf + (»?,rx+ Vy + "zf  z)p.+ (rxrx+ ryr, + rzrz)Pf 

=  -pUCr^u^+ryVf+r.Wf^pVCr^u.+ryV^+r.w,)  (6.8) 

For  simplicity,  the  equation  is  written  as 

aJ>j  +  a2P^  +  a3P^  =  R  ^ /^^^  "     ^     '  (6.9) 

with 

"1  =  (^xrx+Cyry+^zrJ  "    \  . 

"2  =  fex+'J/y+'JirJ  • 
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°3  =  (rxfx+ryr,+fzrz) 

R  =  -pU(5:^u5+ryVj+rzWe)-pV(f,u^+ryV^+r,w^)  (6.10) 

Note  that  the  left-hand  side  of  Eq.(6.8)  equals  vp«Vf  which  is  the  normal 
pressure  gradient  without  normalization. 

The  normal  momentum  equation,  Eq.(6.9),  is  solved  at  the  solid  surface  by 
using  ADI  decomposition  in  the  C  and  ij  directions  and  second-order  forward 
differences  in  the  f  direction.  A  second-order  forward  difference  is  first  applied 
to  Pj.,  i.e., 

Pf  =  (-3Pj,M  +  '*Pjw-Pj,M)/2  +  O(Ar')  (6.11) 
The  normal  momentum  equation,  Eq.(6.9),  can  then  be  written  in  semi-discritized 
form  as 

2ti        2c(  2R 
Pik,r  3;^  Pr  17-  P.  =  ('^Pjw-Pjw-  —  )/3  (6.12) 

If  the  ^-  and  r/-  derivatives  are  central  differenced,  then  the  finite-difference 

equation  of  Eq.(6.12)  expressed  in  operator  form  is 

2a,            2a,                             2R  . 
Pj.k.r  3—«fPj,k,ri— ^.Pjjci  =  (4Pj,M-Pjjc,3  )/3 

or 

^^'^'^'^  '"^Pj'^i  =  (4p^,K,2-Pj,M-  ^  )/3  (6.13) 
By  applying  ADI-de composition  to  the  left-hand  side  of  Eq.(6.13),  one  has 

^        1?-  '"^PJ'M  =  (4Pjw-Pjw-  ^  )/3  (6.14) 

which  is  second-order  accurate  in  each  spatial  coordinate.  The  pressure  at  the 
solid  surface  can  then  be  solved  by  the  following  two  one-dimension  inversions 
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(■|^Vl)P*  =  (4pj^-Pi«-|^)/3 

(^Vl)Pj.M  =  P*  (6.15) 

The  left-hand  side  of  Eq.(6.15)  forms  systems  of  tridiagonal  matrices  which  are 
solved  by  the  Thomas  algorithm. 

The  density  at  the  solid  surface  is  obtained  via  the  ideal  gas  equation  of 
state  when  the  pressure  and  temperature  are  known.  For  an  isothermal  wall  at 
constant  temperature  T^  the  density  is  computed  from  the  dimensionless 
equation  of  state 

P  =  ypTJT^  (6.16) 
For  an  adiabatic  wall,  the  stagnation  enthalpy      is  held  constant.  By  using  the 
equation  for  enthalpy 

H„  =  E+p/(p,p)  =  E„+pJp„  (6.17) 
and  the  computed  velocities  and  pressure,  a  value  of  density  is  obtained  at  the 
solid  surface,  i.e., 

_  yppa 

'  ~   (7-l)[(E.+p,)-,.(u^+v^+w^)/2]  (6-18) 
Viscous  flow.  For  viscous  flows,  the  no-slip  condition  is  enforced  by 
setting  all  the  velocity  components  equal  to  zero  on  the  impermeable  wall,  i.e., 

UjW  =  Yiw  =        =  0  (6.19) 
The  total  energy  is  again  computed  from  Eq.(2.7)  in  which  the  pressure  and 
density  are  obtained  as  follows. 

Since  the  viscous  terms  have  a  negligible  effect  on  the  surface  pressure, 
the  pressure  distribution  for  viscous  flows  is  obtained  by  the  same  procedure  as 
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that  for  inviscid  flow.  It  is  worth  noting  that  for  viscous  flows,  U  and  V  are  both 
zero  on  the  solid  wall;  and  hence,  the  normal  momentum  equation  being  solved 
becomes 

Uxfx+^yfy+e.rJPe+('7xrx+'7yry+'?,fJp,+(r,f«+ryry+fxrz)Pf  =  o  (6.20) 

In  other  words,  based  on  the  inviscid  normal  momentum  equation,  Eq.(6.9),  the 
solid  boundary  condition  for  the  viscous  flow  is 

an  |vr| 

The  density  for  an  isothermal  wall  at  constant  temperature      is  computed  from 
the  same  dimensionless  equation  of  state,  Eq.(6.17).  However,  for  an  adiabatic 
wall,  extrapolation  is  used  for  viscous  flow  problems,  e.g., 

PiM  =  f>w  (6.21) 
since  both  pressure  and  temperature  gradients  are  zero  at  the  wall,  and 

dp 

  =0    at  a  wall 

an 

Slow  start.  In  viscous  flow  problems,  the  flow  field  is  usually  started  from 
an  initial  condition  of  uniform  free-stream  flow.  To  avoid  the  sudden  jump  on 
the  solid  surface,  the  no-slip  condition  is  turned  on  slowly  over  the  first  several 
time  steps.  In  the  code,  the  boundary  conditions  are  smoothed  over  the  first  30 
time  steps  by  setting  the  dependent  variables  q  at  the  sohd  wall  as 

q  =  <^  +  (l-«)q.  (6.22) 
with 

u>  =  (10-155 +  6^')^^  9  =  (NC-l)/30  (6.22a) 
where  NC  is  the  number  of  time  steps  and  q^^  is  the  correct  boundary  conditions 
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discuss  above  [28].  Equation  (6.22a)  is  a  special  function  designed  for  w  such 
that  w  increases  very  slowly  over  the  first  few  time  steps. 

6.2.2  Far-Field  Condition 

,      In  order  to  numerically  solve  a  problem  that  is  posed  in  an  exterior  region, 
the  infinite  domain  is  commonly  replaced  by  a  finite  domain  with  artificial 
bound£uies.  Stretched  grids  are  usually  used  to  place  those  boundaries  far  away 
from  the  object  such  that  no  spurious  reflected  waves  can  be  generated.  Hence, 
the  free-stream  conditions  can  be  specified  as  the  far-field  boundary  conditions 
[30].  •  3    . '  .  - 

The  subroutine  FSBC  is  used  to  impose  the  far-field  boundary  conditions 
at  1=  Imax  by  setting  all  the  dependent  variables  q  at  those  boundaries  equal  to 
free-stream  condition,  i.e.,      .  ..  '  '-'l  i' 

q  =  q«  '  '         '  "         :      '  (6.23) 

6.2.3  Inflow/Outflow  Condition  ;  \ 

The  inflow  and  outflow  boundaries  may  be  viewed  as  special  cases  of  far- 
field  boundaries,  and  hence,  the  non-reflecting  property  is  also  required.  Since 
the  incoming  and  outgoing  directions  are  known  at  those  boundaries,  it  is  easier 
to  specify  the  boundary  conditions  based  on  the  wave  theory  [30].  When  the 
flows  at  the  inflow  boundaries  are  supersonic,  the  flow  conditions  are  not  affected 
by  the  downstream  flow  characteristics  and  hence  can  be  set  equal  to  the  free- 
stream  conditions,  i.e., 

q  =  q. 
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For  supersonic  flow  at  outflow  boundaries,  extrapolation  can  be  employed  since 
no  influence  should  be  felt  upstream  in  the  region  of  interest,  e.g.,  at  j  =  jmax 

4jniax,k44  ~  ^jinax-l,k44 

(6.24) 

For  subsonic  flows  the  situation  is  more  complicated.  More  sophicicated 
boundary  conditions  may  be  needed  if  strong  reflection  waves  are  expected  at 
the  artificial  boundaries.  Here,  only  the  techniques  used  in  the  original  ARC3D 
are  described.  In  ARC3D,  the  free-stream  condition,  Eq.(6.23),  is  employed  at 
inflow  boundaries  while  at  the  outflow  boundaries,  the  velocity  components  u,  v, 
w  and  the  density  p  are  determined  by  extrapolation  and  the  total  energy  E  is 
obtained  from  Eq.(2.7)  by  assuming  the  pressure  is  equal  to  the  free-stream 
pressure,  e.g.,  at  j  =  jmax 

^jinax,k4  ~  Ujmax-lJtJ 
X)max,k,l  ~  yjmax-l,k4 

^jmax,M  ~  ^jmax-l,k4  (6.25) 

E  =  pV(7-1)  +  ^p(u2+v2+w2) 

6.2.4  Symmetric  Boundary  Condition 

The  symmetry  boundary  is  commonly  adopted  to  simplify  flow  problems  if 
symmetric  conditions  can  be  assumed.  In  the  computer  code,  the  »?  =  0  plane  is 
assumed  to  be  a  symmetry  plane  which  is  also  a  constant  z  plane  in  physical 
space  as  shown  in  Figure  6.5a.  On  this  plane,  the  V-velocity  and  the  first 
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(a)  Symmetric  boundary 


Synmetry  Plane 
Y  z-0 

Mirror  image  plane 

k-32 

2-0 

A.  

1 

Mirror  imoge  plane 

(b)  Symmetric  boundary  with  mirror  image  plane 
Figure  6.5  Symmetric  boundary  condition 
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derivatives  of  the  other  flow  properties  with  respect  to  n  are  set  to  zero  to 
enforce  symmetry,  i.e.,     qj      =  0 

%W  =  (Hmw  -  qj,k+ui)/3  for  i  ^  4  (6.26) 
Recall  that  qj^y^  is  pw  at  grid  point  and  V=w  on  this  symmetry  plane. 

Another  way  of  dealing  with  symmetry  boundaries  is  to  add  one  more  plane  that 
is  the  mirror  image  of  the  plane  next  to  the  symmetry  plane.  For  instance,  if  k  = 
2  is  the  symmetry  plane,  then  k  =  1  is  the  added  boundary  plane  which  is  the 
mirror  image  of  k  =  3  plane.  As  shown  in  Figure  6.5b,  if  the  k =2  plane  is  also  a 
constant  z  plane  then  the  flow  conditions  at  k  =  1  plane  can  then  be  specified  as 

%W  =  q^Ai       for  i  ^  4  (6.27) 
to  enforce  symmetry.  Note  that  q^y^  is  (pw)j^y  and  equals  -(pw)^^  since  k  =  1 
plane  is  a  constant  z  plane. 

6.3  Interface  Boundary  Condition 

The  success  of  the  multi-block  grid  method  relies  on  correct 
communication  among  blocks.  Inappropriate  implementation  of  interface 
boundary  conditions  may  result  in  numerical  instability  or  reduced  accuracy. 
There  have  been  various  techniques  developed  for  updating  the  interface 
boundary  conditions  with  different  degree  of  sophicication  [52-53].  A  simple  and 
efficient  interpolation  method  has  been  successfully  tested  on  a  projectile  flow 
problem.  In  this  approach,  the  information  at  the  interface  boundaries  is 
obtained  by  distance  weighted  interpolation  from  the  neighboring  blocks.  For 
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instance,  the  unknowns  q  at  an  interface  boundary  point  can  be  estimated  from 
the  known  quantities     at  an  interior  point  at  block  1  and  q2  at  block  2  as 
follows: 

q  =  (qiSi  +  q2Si)/(si  +  &z)  (6.28) 

where  s^  and  Sj  are  the  distances  of  the  interior  points  to  the  boundary  point.  In 
other  words,  the  gradient  of  q  in  the  local  one-dimension  direction  s  is  assmned 
constant,  i.e.. 

It  is  noted  that  the  conservative  properties  are  not  conserved  in  this  method  and 
errors  may  be  introduced  if  the  solution  variations  are  large  near  the  interfaces. 
However,  the  method  is  quite  efficient  and  easy  to  implement  and  has  been 
tested  on  a  projectile  flow  problem  with  satisfactory  results.  Numerical 
experiments  showed  that  only  negligible  errors  were  introduced  if  the  grid  was 
smooth  and  orthogonal  near  the  interfaces. 

The  interfaces  of  the  six-block  flow  domain  are  identified  by  the  block 
diagrams  shovm  in  Figure  6.6.  Here,  some  special  lines  that  require  special 
treatment  are  discussed.  As  shown  in  Figure  6.6a,  the  line  j=28,  k=13  of  Block  1 
is  a  common  line  among  Block  1,  2,  3,  4  and  5,  the  boundary  conditions  at  this 
line  are  obtained  by  interpolating  over  all  these  blocks,  i.e.,  for  1  =  1,  60 

qi28434  =  q2i,i24  =  q3i,ij  =  q44o,i4  =  ^^iw 

=  (Slj27q32,u  +  S3j2ql27,i34  +  ^^i^<l540W  +  ^^^^ihu^)/ 
(Slj27  +  S3j2  +  S2,j2  +  s5u2)  (6.30) 
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(a)  Block  1  &  2,5  and  Block  6  &  2,5 
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(b)  Block  2&3,  3&4and4&5 


Figure  6.6  Interface  block  diagram 
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Here,  ql28,u4  stands  for  the  unknown  q  at  grid  point  (28,13,1)  of  block  1  and  slj27 
stands  for  the  distance  between  grid  point  (27,13,1)  and  (28,13,1)  in  block  1,  and 
likewise  the  other  notation  can  be  defined.  A  similar  formula  is  applied  to  the 
line  j  =  1,  k=  13  of  Block  6  which  is  a  conmion  line  of  Block  2,  3,  4,  5  and  6. 
As  can  be  seen  from  the  grid  network,  the  three  surfaces  j  =  l,  k=35  and  j=40  of 
Block  3,  which  are  also  j=  40,  k=35,  j  =  l  surfaces  of  Block  4,  form  almost  a  flat 
plane.  The  boundary  conditions  at  the  lines  between  the  surfaces,  i.e.,  j  =  l,  k=35 
and  j=40,  k=35  of  Block  3,  are  interpolated  not  only  from  block  3  and  4  but  also 
from  the  neighboring  surfaces,  i.e.,  for  1  =  1,  60 

<l3i^54  =  QVw  =  (s3k34q32^54  +  S3j2q3i^^/(s3k34  +  s3j2)  * 

q34o,354  =  =  (s4k34q3„s4  +  S4j2q3i^^)/(s4k34  +  s4j2)         '  ;    .  (6.31) 

6.4  Convergence  Criterion  :-}:r-.  - 

The  numerical  algorithm  employed  in  the  code  is  a  time-marching  scheme. 
The  steady-state  solution  can  be  obtained  as  a  time  asymptotic  solution  of  the 
unsteady  Navier-Stokes  equations.  If  one  rewrites  the  numerical  scheme, 
Eq.(3.39),  as  follows 

[I + xes    + ]"[I + Ae5,B + D  J° 

[I+Ae(6^C-Re-^5f(J-'MJ))  +  D^]"Aq"-DE  =  RHS 

RHS  =  -A(5jE+5/+5^G-Re-^Sj.S)" 
the  vector  RHS  will  approach  zero  if  the  computations  approach  a  steady-state 
solution.  Hence,  the  rate  of  convergence  can  be  indicated  by  the  averaged  12 
residual  computed  as  follows 
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^  "'^'^  '  A,[amJlKS')(lmax-l)l"  <"2' 
Here,  jmax,  kmax,  and  Imax  are  the  number  of  grid  points  in  each  curvilinear 
directions. 


6.5  Pressure  and  shear  stress  coefficient 
The  pressure  and  shear  stress  coefficients  which  are  defined  as 

_   _        r  _  r'/Re 

are  computed  as  follows: 

(a)  pressure  coefficient 

p*  is  first  calculated  from  Eq.(2.9)  and  p*,  =  l/y,  Cp  is  then  readily 
computed  from  Eq.(6.33). 

(b)  shear  stress  coefficient 

The  shear  stress  tensor  in  Cartesian  coordinate  system  can  be  obtained 
firom  Eq.(2.2b).  Since  only  the  shear  stresses  on  the  solid  surface  are  of  interest, 
the  shear  stress  tensor  is  then  expressed  in  the  curvilinear  coordinate  system,  i.e., 

[r(^,n,^)]  =  [R][r][R]-»  (634) 
where  [R]  is  the  transformation  matrix  that  can  be  expressed  in  terms  of  the 
normaUzed  metrics  coefficients.  The  shear  stress  components  t^^  and  t^^  are 

r,,  =  (Ae,+Bey+CO/(^,'+c/+e/)^  (6.35) 

r,,  =  (Ar,,  +  Bv,+  CrjJ/(r,,'  +  v/+V.Y  (6,35a) 
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where 

B  =  (r^^  +  r^y^  +  r^,)/(x,2+y^'+z^')^ 

C  =  ('■,^f  +  ryj,  +  r^,)/(x/+y,^+z,^)^  (6.35c) 
The  shear  stress  coefficients  are  then  computed  from  Eq.(6.35a),  i.e., 
Cffj  =  2r5.j/ReM„ 

Cf,^  =  Ir'yReM^        •  (6.36) 
The  drag  and  Uft  forces  can  be  obtained  by  integration  of  the  surface  pressure 
and  shear  stress  distributions  over  the  body  surface. 


CHAPTER  VII 
RESULTS  AND  DISCUSSION 


The  proposed  multi-block  computational  method  is  investigated  for  the 
simulation  of  transonic  turbulent  flows  of  Mach  0.8  past  a  wing-fuselage 
configuration  by  using  the  six-block  grid  network.  Since  the  only  available 
experimental  data  are  the  surface  pressure  coefficient  at  eight  wing  cross-sections, 
the  discussion  will  concentrate  on  the  results  computed  at  these  locations.  These 
eight  wing  cross-sections,  expressed  in  normalized  distances  from  the  symmetry 
plane,  are,  respectively,  0.15,  0.225,  0.305,  0.405,  0.5,  0.7,  0.85,  and  0.95,  in  which 
the  distance  from  the  symmetry  plane  to  wing  tip  is  one. 

The  multi-block  thin-layer  Navier-Stokes  code  (WF6B)  provided  the 
options  to  compute  the  Jacobian  and  metric  coefficient  by  either  finite-difference 
approach  or  volume/cell-face  area  approach.  The  use  of  finite-difference 
approach  for  computing  Jacobians  of  transformation  results  in  negative  values  at 
the  certain  regions  near  wing  tip,  the  leading  edge,  and  trailing  edge.  This  seems 
to  indicate  that  the  grid  network  is  quite  skew  near  these  regions.  Hence,  in 
order  to  accurately  compute  the  geometric  quantities,  it  is  necessary  to  use  the 
volume  approach  for  Jacobian  calculation.  The  metric  coefficients  can  be 
computed  by  either  approach.  It  is  known  that  the  use  of  the  cell-face  area 
approach  for  computing  metric  coefficient  is  more  consistent  with  the  volume 
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approach  for  Jacobian  calculation;  however,  the  computed  pressure  coefficients 
are  found  to  be  insensitive  to  either  method. 

7.1  The  Result  from  Base  Grid 

The  computational  method  was  first  investigated  for  a  transonic  mrbulent 
flow  of  Mach  0.8  past  the  wing-fuselage  configuration  at  zero  angle  of  attack  with 
the  use  of  the  coarse  base  grid.  The  grid  sizes  of  block  1  to  block  6  are 
28x25x60,  40x13x60,  40x35x60,  40x35x60,  40x13x60,  and  26x25x60,  respectively.  A 
converged  solution  was  obtained  in  2500  time  steps.  Figure  7.1  shows  the 
comparison  of  the  computed  Cp-distribution  and  the  experimental  data  at  four 
wing  cross-sections.  The  result  shows  that  the  calculated  wing  surface  pressure 
coefficient  agrees  well  qualitatively  with  the  experimental  data  but  the 
quantitative  agreement  is  not  satisfactory  at  all  at  certain  locations.  One  can  " 
clearly  observe  the  large  discrepancy  in  the  leading  edge  area  and  around  the 
shock  wave  and  high  pressure  gradient  regions.  In  order  to  investigate  the  effect 
of  the  leading  edge  grid  resolution  upon  solution  accuracy,  the  base  grid  was 
modified  by  adding  10  more  grid  points  in  the  regions  close  to  the  leading  edge 
of  the  wing.  This  revised  grid  system  was  then  used  for  a  computation  that 
restarted  from  the  afready  obtained  solution.  The  computed  Cp-distribution 
obtained  in  500  time  steps  is  shown  in  Figure  7.2,  and  does  not  show  any 
improvement.  This  seems  to  indicate  that  around  the  leading  edge,  the  grid 
spacing  is  not  the  major  cause  of  the  discrepancy.  It  is  suspected  that  the  highly 
skewed  grid  in  that  region  as  well  as  the  turbulence  model  may  be  responsible  for 
the  disagreement  between  the  prediction  and  the  measurement. 


Figure  7.1  Cp-plot  for  a=0°  based  on  the  base  grid 
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Figure  7.2  Cp-plot  based  on  leading  edge  refined  grid 
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In  order  to  gain  further  insight  and  understanding  on  the  solution 
algorithm,  computation  was  made  for  flow  at  angle  of  attack  a =-3°  by  using  the 
base  grid.  The  solution  obtained  from  the  case  of  zero  angle  of  attack  was  used 
as  an  initial  guess  for  restarting  the  computation.  The  Cp-distribution  obtained 
after  400  time  steps  is  shown  in  Figure  7.3.  As  shown  in  the  figure,  the  Cp- 
distribution  on  the  upper  wing  surface  agrees  well  with  the  measured  data; 
however,  the  agreement  is  not  acceptable  on  the  leeward  side  of  the  wing.  This 
seems  to  suggest  that  the  grid  resolution  is  not  fine  enough  to  resolve  the  flow 
phenomena  involving  shock  waves  and/or  large  solution  variations.  Hence,  it  is 
demanded  for  a  grid  refinement  study. 

7.2  The  Result  from  Refined  Grid 

Limited  by  computer  resources  and  accessibility,  only  the  grid  resolution  in 
the  wing  chordwise  direction  was  improved.  A  refined  grid  was  generated  by 
using  natural  cubic  spline  to  double  the  number  of  grid  points  in  the  wing 
chordwise  direction.  With  the  use  of  the  refined  grid,  flow  at  zero  angle  of  attack 
was  again  computed.  A  converged  solution  was  obtained  in  500  time  steps  after 
restarting  from  the  result  of  the  base  grid.  Figure  7.4  presents  the  calculated  Cp- 
distribution  and  the  experimental  data  at  eight  different  cross-section  along  the 
wing  span.  As  can  be  seen,  the  computed  surface  pressure  coefficient  indeed 
agrees  more  closely  to  the  measured  data  now,  particularly  in  the  shock  and  high 
gradient  regions.  These  results  clearly  indicate  that  the  accuracy  of  numerical 
sunulation  can  be  further  improved  if  a  further  improved  grid  resolution  or  a 
solution-adaptive  grid  is  used  in  the  simulation.  Figure  7.5  shows  the  pressure 


Figure  7.4  Cp-plot  for  0=0°  based  on  refined  grid 


Figure  7.4  (continued) 
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(b)  Top  and  side  views;  Back  fuselage 
Figure  7.5  Cp-contour  for  a=0" 
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(c)  Top  and  side  views;  Block  5  fuselage 


(d)  Top  and  side  views;  Block  2  fuselage 


Figure  7.5  (continued) 
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0 


(e)  Upper  wing  surface 


(f)  Lower  wing  surface 
Figure  7.5  (continued) 
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coefficient  contour  on  the  surface  of  the  wing-fuselage  configuration  with  an 
interval  of  ACp=0.05.  A  Mach  contour  plot  on  a  surface  1=20,  a  short  distance 
above  the  wing  surface,  is  given  in  Figure  7.6,  which  provides  the  information  of 
the  location,  extent  and  strength  of  shock  wave.  The  distributions  of  calculated 
shear  stress  components,  t^^  and  t^^  on  the  wing  surface,  as  given  in  Figures  7.7- 
7.9,  show  that  there  is  no  reversed  flow  in  the  streamwise  direction  on  the  wing 
and  that  the  fluid  particles  sUghtly  above  the  wing  surface  are  moving  away  from 
the  wing  root  toward  the  wing  tip  except  that  some  of  those  near  the  under 
surface  (block  3)  trailing  edge  move  toward  the  wing  root  as  they  are  swept 
downstream.  It  was  observed  that  the  flow  characteristics  represented  by  the 
shear  stress  distribution  are  in  agreement  with  those  exhibited  by  the  surface 
pressure  distribution.  For  instance,  the  variation  of  r^^  shown  in  Figure  7.7 
clearly  indicates  a  sudden  decrease  in  fluid  velocity  across  the  shock. 

To  further  investigate  the  method,  the  refined  grid  was  then  used  for  the 
simulation  of  flows  at  4"  angle  of  attack.  The  result  obtained  from  the  base  grid 
was  again  used  as  the  initial  guess  for  the  computation.  Figure  7.10  presents  the 
Cp-distribution  obtained  in  500  time  steps.  The  computed  Cp  on  the  windward 
side  is  in  good  agreement  with  measure  data;  however,  the  agreement  on  the 
leeward  is  not  favorable.  The  predicted  shock  location  on  the  upper  surface  is 
downstream  of  the  measured  location.  It  is  also  noticed  that  the  computed  Cp 
has  shown  oscillations  near  the  leading  and  trailing  edges  of  the  wing.  It  is 
argued  that  the  oscillations  could  be  caused  by  the  improper  use  of  the  cubic 
spline  technique  in  refining  the  grids,  since  they  appear  only  in  the  results 
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(a)  Upper  wing,  block  4 


(b)  Lower  wing,  block  3 


Figure  7.6  Mach-contour  at  1=20 


Figure  7.7  r^^-plot  for  0=0"  based  on  refined  grid 
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(b)  Lower  wing  surface 
Figure  7.8  Surface  shear  stress  component  r 
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(a)  Upper  wing  surface 


(b)  Lower  wing  surface 
Figure  7.9  Surface  shear  stress  component  r 


■1 

I 

t 


127 


Figure  7.10  Cp-plot  for  0=4" 
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obtained  from  the  refined  grid.  The  use  of  cubic  spline  may  result  in  more 
skewed  grid  from  the  original  skew  base  grid.  The  calculated  shear  stress 
distribution  is  shown  in  Figure  7.11,  indicating  that  there  is  shock  induced 
separation. 

Finally,  some  information  regarding  the  CPU  time  requirement  for  the 
problems  studied  is  given  to  help  illustrate  the  scope  as  well  as  the  complexity  of 
the  present  methodology.  All  the  computations  were  performed  at  the 
supercomputer,  Cray-2  of  NAS  system  at  NASA  Ames  Research  Center.  The 
computer  code  WF6B  has  the  efficiency  of  about  7.7x10'^  second  per  grid  point 
per  time  step.  The  base  grid  has  totally  311,400  grid  points,  and  requires  about 
24  second  for  one  time  step.  The  result  for  zero  angle  of  attack  were  computed 
for  2,500  time  steps,  that  is,  about  17  hours  of  CPU  time.  The  refined  grid  has 
536,040  grid  points  and  requires  41  second  for  one  time  step.  The  result 
obtained  was  computed  for  500  time  steps  after  restarting  from  the  result  of  the 
base  grid.  The  CPU  time  is  about  5.5  hours  for  500  time  steps. 


Figure  7.11  Shear  stress  coefficient  plot  for  0=4" 


CHAPTER  Vm 
CONCLUDING  REMARKS 

A  multi-block  computational  method  is  developed  for  three  dimensional 
high  speed  turbulent  flows  over  complex  configurations.  In  this  method,  the  flow 
field  is  divided  into  several  contiguous  blocks  such  that  each  block  is  governed  by 
the  same  set  of  thin-layer  Navier-Stokes  equations.  The  solution  algorithm  with 
distance  weighted  interpolation  for  updating  the  interface  boundary  condition  is 
rather  simple  and  straightforward;  however,  special  measures  in  the  grid  topology 
and  grid  generation  are  required  before  a  satisfactory  prediction  can  be  made.  A 
flow  simulation  package  which  includes  grid  generation,  flow  solver,  and  plotting 
program,  is  developed  and  apphed  to  a  transonic  turbulent  flow  of  Mach  0.8  past 
a  wing-fuselage  configuration.  The  flow  solver  is  good  for  high  Reynolds  number 
transonic  turbulent  flow  problems  with  modest  flow  separations.  A  special 
designed  six-block  grid  topology  is  used  for  the  development  of  the  method.  The 
flow  field  is  properly  divided  into  six  blocks  that  are  solved  by  the  same  unsteady 
TVD  thin-layer  Navier-Stokes  code.  A  coarse  base  grid  and  a  refined  grid  are 
generated  and  used  for  flow  simulations. 

Numerical  experiments  performed  in  this  study  have  shown  that  special 
measures  and  care  are  required  in  generating  good  grid  system  involving  excessive 
distortion  between  the  physical  and  computational  domain;  however,  the 
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computed  shock  location  and  pressure  distribution  can  be  in  good  agreement  with 
the  experiment  data  if  the  block  grids  used  are  adaptive  to  important  solution 
characteristics.  The  proposed  method  has  demonstrated  its  flexibility,  efficiency, 
and  capability  in  the  simtilation  of  complex  geometry  flow  problems.  In  the 
demand  of  the  CFD  methods  for  simulating  turbulent  flows  over  complex 
geometry,  the  proposed  method  is  indeed  a  very  promising  approach  to  be  further 
developed. 

As  experienced  in  this  study,  a  good  grid  network  is  crucial  to  achieve 
accurate  results.  To  develop  a  cost-effective  computational  method  for  3-D 
complex  configiiration  aerodynamics,  a  technique  to  generate  solution-adaptive 
grid  network  is  desired.  However,  a  good  grid  that  adapts  to  the  important  flow 
characteristics  of  the  present  complex  flows  are  not  known  a  priori;  hence,  it  is 
very  desirable  to  incorporate  a  self-adaptive  grid  generation  technique  that  can 
give  a  proper  grid  distribution  for  yielding  accurate  solutions  with  minimum 
number  of  grid  points  to  minimize  the  computing  cost. 


APPENDIX  A 
JACOBIAN  MATRICES 

The  Jacobian  matrices  A,  B  or  C  equals 

k.  kx  k,  0  ' 

k,^-u<?  k,+  (9-(7-2)k,u  kyU-(7-l)k,v  k,u-(7-l)k,w  k,(7-l) 

ky^-v(?  k,v-(7-l)kyU  k,+  (?-(7-2)kyVk,v-(7-l)kyW  ky(7-l) 

k^^-wfl  k,w-(7-l)k^u  kyW-(7-l)k,v  k,+  <?-(7-2)k^w  k,(7-l) 

.  e{4,-^)  kj>-(cf-\)u6  ky^-(y-l)ve    kji,-(-r-l)we     k^+js  . 

where 

i>  =  7(e/p)-<^ 

<f>  =  (7-l)(u^+v^+w^)/2 

e  =     +  ky  + 
with  k  =  ^,  T7  or  r  for  A,  B  or  C,  respectively. 
The  viscous  flux  Jacobian  M  is 


0 

0            0            0  0 

«l«f(Vp)    a2«f(l/p)    «3«f(Vp)  0 

"2«f(l/p)     «4«f(Vp)  0 

«35f(l/p)    «5«f(l/p)    c^eScii/p)  0 

.  ^51 

mjj                        a„5f(l/p)  , 

where 

ni2i 

=       ^(-u/p)  +  a^S  j.(-v/p)  +  ajSf  (-w/p) 
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=  a2«f(-u//')  +  a45f(-v/p)  +  a5«f(-w/p) 

ni4i  =  a3«f(-u/p)  +  Q5Sf(-v//')  +  a6«f(-w/p) 
nisi  =  ai6f(-uVp)+a25f(-2uv/p)+«35f(-2uw/p) 
+  a^S  f  (V/p  )  +  aj6  ^(  V/p)  +  aj*  ^  (-2vw/p) 

°l52  =  -ni2i-ao*f(u/p) 

=  -m4i-ao5f(w/p) 

ao  =  7MPr-'(r/+ry'+rz') 
«i  =  4(4/3)r/+fy'+r/] 
"2  =  W3)r.ry 

«4  =  M[f,'+(4/3)r/+f,'] 

"5  =  (M/3)ryrz 

«6  =  /^[r,'+f/+(4/3)r,'] 


APPENDIX  B  : 
EIGEN-FUNCnONS  ; 

The  eigenvalues  of  A,  B  and  C  are 

{a/,a/,a/,a,^a/}  =  {U,U,\J,lJ+ci^,'+i/+^,y,V-c(i,'+i/+Cf} 
{a,^a,^a/,a,^a/}  =  {V,V,V,V^-c(,,^+,/+,,^)^V-c(,,2+,/+,,^)^}  - 
{a/,a,^a,^a,^a,^}  =  {W,W,W,W+c(f,2+fy2+r/)^W-c(r,'+r/+rz')^} 

respectively,  where  c  =  (yp/p)^  is  the  local  speed  of  sound  and  the  contravariant 

velocity  components  U,  V,  W  are 

U  =  C,+l,u  +  CyV+^^  *  '  ^    ^     'V  . 

V  =  i7t  +  >7xU  +  »7yV+i7,W  ■  "t'? 

w  =  c,+rxU+r,v+f^  •  ■  ^rv.      ■>  ■ 

The  similarity  transformation  for  the  Jacobian  matrices  A,  B,  C  are 
A  =  R^,R,  i  ,  B  =  R^A^R^  i  ,  C  =  R,A,R,-^ 
where  the  a's  are  the  diagonal  matrices 

A,  =  D[u,u,u,u+c(^/+e/+e,2)^u-cU/+?/+?,2)«]  :.,  y  v 

A,  =  D[V,V,V,V+c(„^+„/+„/)^V-c(„/+„/+„,2)V4j 

A,  =  D[w,w,w,w+c(^/+f/+^,')^w-c(f/+r/+f/)^] 

If  one  defines  the  notation  as  follows  and  k  =  ^,  »;  or  f  * 
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e  =  ICjU  +  ICyV  +  lCjW 

and 

a  =  J2p/2c,     ^  =  /Z/2pC,    7  =  7-1 

0  =  (7-l)(u^+v2+w2)/2 
then  = 


a 

a 

a(u  +  lc,c) 

q(u-ICjC) 

lc,v-Jc,p 

a(v  +  lCyC) 

a(v-ICyC) 

a(w  +  I^c) 

a(w-I^c) 

p(k,w-I^u) 

p(RyU-Ic^v) 

c[(<l>  +  <^)/j 

+  ce] 

a[(^  +  c')/7 

-ce] 

and  R^-*  = 

^(i-^/<):  ^S'^u/c^-ic./p  ^7v/c^  i^7w/c^+E,/p 

(lc,w-R,u)/p 

R,(lWc^)-  ^yu/<^+yp  1c,7v/c^-Vp  1^7w/c^  -1^7/c' 

^(^-Ce)  ^[Ic,C-7u]  ^[lSC-7V]  ^[1^C-7W]  /37 

^P(<f>  +  Ce)  -p\^C  +  yu]  -^[EyC  +  7V]  -/3[I^c  +  7w]  ^7 


APPENDIX  C 
CLUSTERING  FUNCTIONS 


The  algebraic  grid  generation  requires  some  sort  of  clustering  functions  or 
stretching  functions  to  distribute  the  grid  points  so  that  more  grid  points  are 
clustered  where  the  solution  varies  rapidly,  e.g.  near  the  soUd  surface.  Four 
clustering  functions  have  been  implemented  in  GG3D  and  are  discussed  as  follows: 
1.  CUBICLU  -  Cubic  Function  Distribution 

For  a  space  curve  in  the  1-direction,  for  instance,  the  normalized  arc  length 
s  is  defined,  in  terms  of  arc  length  S(l),  as 
S(l)-S(l) 

^®  =  S(lma.)-S(l)  (^1' 
Let  the  specified  first  spacing  and  the  last  spacing  be  denoted  by  asj  and  AS2, 
respectively.  The  cubic  function  distribution  is  obtained  as  follows. 
First,  define 
h  =  L/(hnax-l) 
Then  a,  b,  c  are  computed  as 
a  =  [asi  +  ASj  -2h] 

b  =  [ASi  -  h  -  c(h^-h)]/(h^-h)  (C2) 
c  =  1  -  b  -  c 

The  normalized  cubic  function  distribution,  v(l),  is  then  given  by 

v(l)  =  a(l-l)  +b(l-l)^  +  c(l-l)'  (C3) 
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2.  TANHCLU  -  Hyperbolic  Tangent  Distribution 

Using  the  same  notation  as  above,  the  hyperboUc  tzingent  distribution  is 
constructed  as  follows. 

First,  compute  A  and  B 

A  =  Vasj/Vasi 


B  =  (lmax-l)yAS2ASi  (C.4) 
Then  the  following  nonlinear  equation  is  solved  for  6  by  the  Newton-Raphson 
method. 

sinh(5)/5  =  1/B  (C.5) 
The  normalized  hyperbolic  tangent  distribution,  v(l),  is  then  given  by 

=  Arfi-A)u(i) 

where 

urn  =  -a+  tanh[5{(l-l)/(hnax-l)-0.5}] 

2  ^  tanh(fi/2)  ^  ^^''^ 

3.  SINHCLU  -  Hvperbolic  Sine  Distrhutinn 

If  only  the  first  spacing  asi  is  specified,  the  hyperbolic  sine  distribution  is 
calculated  as  follows. 

First,  B  is  computed 
B  =  (hnax-l)ASi  (C.8) 
and  6  is  determined  by  the  Newton-Raphson  method  as  the  solution  of 

sinh(5)/5  =  1/B  ,  ■ 

The  normalized  hyperbolic  sine  distribution,  v(L),  is  given  by  ? 
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smh[s{(l-l)/(lmax-l)}]  ^ 
sinh(6) 

4.  EXPOCLU  -  Exponential  Distribution  '  ^■ 

If  only  the  first  spacing,  asj,  is  specified,  the  exponential  distribution  is 
calculated  as  follows. 

First,  5  is  obtained  by  the  Newton-Raphson  method  fi"om  the  following 
nonlinear  equation: 

1  =  ASi[(l  +  5)'™"-*-l]/«  '    '':     ■  (CIO) 

Then  the  normalized  exponential  distribution  is  given  by  .  ^  .      i  * 

_  ,  -  ■  i> 

V(l)  =  V(l-l)  +  ASi(l  +  5)"  V     .  (C.11) 

. ,     ^  —        :  "  ^  . 

5.  Summary  - 

The  truncation  error  is  strongly  affected  by  the  point  distribution,  and  studies 
of  the  effect  of  different  distribution  functions  on  the  truncation  error  have  been 
made  in  the  past.  The  exponential  distribution  can  cluster  more  grid  points  near  the 
1=1  surface.  However,  the  variation  of  the  spacing  is  large.  The  cubic  function  is 
good  for  clustering  points  near  both  ends.  The  hyperbolic  sine  function  gives  a 
smoother  distribution  in  the  immediate  vicinity  of  the  1  =  1  surface,  while  the 
hyperboHc  tangent  function  gives  a  smoother  overall  distribution. 
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